Load packages

library(ggplot2)
library(cowplot)
library(lme4)
library(pbkrtest)
library(tidyr)
library(dplyr)
library(optimx)
library(broom)

Data description

Five censuses (15fa, 16sp, 16fa, 17sp and 17fa) were conducted. Information of neigboring trees (area of conspecific, heterospecific and overall species) 20 radius around the focal quadrats are added to the dataset.

Survival analysis

getwd()
## [1] "E:/R/My code/Git/Changbaishan/Changbaishan1"
rm(list=ls())

## Read in data

pestdat.ag_adult <- read.csv("data/pestdat.ag_adult.csv")
pestdat.ag_adult <- pestdat.ag_adult[,-c(1)]
pestdat.ag_adult$pesticide <- factor(pestdat.ag_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pestdat.ag_adult$census <- factor(pestdat.ag_adult$census)
pestdat.ag_adult$exclosure <- factor(pestdat.ag_adult$exclosure)
pestdat.ag_adult$site <- factor(pestdat.ag_adult$site)
summary(pestdat.ag_adult)
##  site        quadrat           sp.       pesticide exclosure growth.form 
##  1:1732   Min.   :101.0   FRMA   : 765   W:1087    0:2269    shrub:1741  
##  2:1620   1st Qu.:202.0   TIAM   : 673   F:1147    1:2063    tree :2591  
##  3: 980   Median :304.0   PHSC   : 327   I:1065                          
##           Mean   :304.9   EUMA   : 302   C:1033                          
##           3rd Qu.:409.0   ACPS   : 296                                   
##           Max.   :510.0   ACBA   : 260                                   
##                           (Other):1709                                   
##   census         deaths            total            survs       
##  16fa:1169   Min.   : 0.0000   Min.   : 1.000   Min.   : 0.000  
##  16sp: 625   1st Qu.: 0.0000   1st Qu.: 1.000   1st Qu.: 1.000  
##  17fa:1406   Median : 0.0000   Median : 2.000   Median : 1.000  
##  17sp:1132   Mean   : 0.5919   Mean   : 2.795   Mean   : 2.203  
##              3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.: 3.000  
##              Max.   :16.0000   Max.   :23.000   Max.   :20.000  
##                                                                 
##   quad.unique    plot             quad.sp         A.sum      
##  1_1_407:  37   1_0:859   1_0_101_ACKO:   4   Min.   : 3990  
##  1_1_503:  36   1_1:873   1_0_101_CEOR:   4   1st Qu.: 4991  
##  2_1_405:  35   2_0:864   1_0_101_EUMA:   4   Median : 5452  
##  1_0_510:  34   2_1:756   1_0_101_FRMA:   4   Mean   : 5631  
##  2_0_510:  33   3_0:546   1_0_102_CEOR:   4   3rd Qu.: 6084  
##  1_1_408:  32   3_1:434   1_0_102_DEGL:   4   Max.   :10411  
##  (Other):4125             (Other)     :4308                  
##      A.con              A.het          con.dens       A.con_curt    
##  Min.   :   0.000   Min.   : 2374   Min.   : 0.00   Min.   : 0.000  
##  1st Qu.:   0.000   1st Qu.: 4487   1st Qu.: 0.00   1st Qu.: 0.000  
##  Median :   4.235   Median : 5182   Median : 1.00   Median : 1.618  
##  Mean   : 427.093   Mean   : 5204   Mean   : 5.41   Mean   : 4.295  
##  3rd Qu.: 658.840   3rd Qu.: 5764   3rd Qu.: 9.00   3rd Qu.: 8.701  
##  Max.   :5308.942   Max.   :10411   Max.   :46.00   Max.   :17.445  
##                                                                     
##    A.het_curt      A.sum_curt    con.dens_curt        dens         
##  Min.   :13.34   Min.   :15.86   Min.   :0.000   Min.   :-0.30365  
##  1st Qu.:16.49   1st Qu.:17.09   1st Qu.:0.000   1st Qu.:-0.30365  
##  Median :17.30   Median :17.60   Median :1.000   Median :-0.13447  
##  Mean   :17.23   Mean   :17.74   Mean   :1.044   Mean   : 0.00000  
##  3rd Qu.:17.93   3rd Qu.:18.26   3rd Qu.:2.080   3rd Qu.: 0.03472  
##  Max.   :21.84   Max.   :21.84   Max.   :3.583   Max.   : 3.41845  
## 
### Intial plots
ggplot(data=pestdat.ag_adult, 
       aes(x=pesticide, y=survs/total, colour=as.factor(exclosure))) + 
  geom_boxplot() +
  facet_wrap(~growth.form)

group_by(pestdat.ag_adult, pesticide, exclosure, growth.form) %>% 
  summarise(surv = mean(survs/total), 
            se.survs=sd(survs/total)/sqrt(length(survs/total))) %>% 
ggplot(aes(x=pesticide, y=surv, ymin=surv-se.survs, ymax=surv+se.survs, 
           colour=as.factor(exclosure))) + 
  labs(x="Pesticide", y="Survival") + 
  geom_errorbar(position=position_dodge(width = .5), width = 0.2) + 
  geom_point(position=position_dodge(width = .5)) +
    facet_wrap(~growth.form)

###------------------------- Initial Survival models -------------------------###
pestdat.ag_adult$site <- factor(pestdat.ag_adult$site)
contrasts(pestdat.ag_adult$site) <- "contr.sum"

## Conspecific and heterospecific area of adult trees
mod.tree_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))


# sum area of adult trees
mod.tree_adult2 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.sum_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

## conspecific density of adult trees
mod.tree_adult3 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + con.dens_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult1, mod.tree_adult2, mod.tree_adult3)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult2: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult2:     A.sum_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult3: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult3:     con.dens_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult1:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                 Df    AIC    BIC  logLik deviance   Chisq Chi Df
## mod.tree_adult2 13 5302.2 5378.4 -2638.1   5276.2               
## mod.tree_adult3 13 5299.0 5375.2 -2636.5   5273.0  3.1694      0
## mod.tree_adult1 14 5269.1 5351.2 -2620.6   5241.1 31.9072      1
##                 Pr(>Chisq)    
## mod.tree_adult2               
## mod.tree_adult3  < 2.2e-16 ***
## mod.tree_adult1  1.617e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## interaction terms
mod.tree_adult4 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure*census +
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))


mod.tree_adult5 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure + census +
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))


mod.tree_adult6 <- glmer(cbind(survs, deaths) ~ site + (pesticide + exclosure)* census +
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult4, mod.tree_adult5, mod.tree_adult6)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult5: cbind(survs, deaths) ~ site + pesticide * exclosure + census + 
## mod.tree_adult5:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult6: cbind(survs, deaths) ~ site + (pesticide + exclosure) * census + 
## mod.tree_adult6:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult4: cbind(survs, deaths) ~ site + pesticide * exclosure * census + 
## mod.tree_adult4:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.tree_adult5 17 5272.3 5371.9 -2619.1   5238.3                         
## mod.tree_adult6 26 5254.2 5406.6 -2601.1   5202.2 36.009      9  3.951e-05
## mod.tree_adult4 38 5268.3 5491.0 -2596.2   5192.3  9.938     12     0.6214
##                    
## mod.tree_adult5    
## mod.tree_adult6 ***
## mod.tree_adult4    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# without neighoring trees
mod.tree_adult7 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
                           (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))


# nested random effects
mod.tree_adult8 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + 
                           A.con_curt + A.het_curt +
                     (1|plot/quadrat) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult1, mod.tree_adult8)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree_adult8:     A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult1:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.tree_adult8 13 5269.0 5345.2 -2621.5   5243.0                         
## mod.tree_adult1 14 5269.1 5351.2 -2620.6   5241.1 1.8897      1     0.1692
mod.tree_adult9 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + 
                           A.con_curt + A.het_curt +
                     (1|plot/quadrat) + (1|plot/census) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
mod.tree_adult10 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +
                           A.con_curt + A.het_curt +
                     (1|plot/pesticide/quadrat) + (1|plot/census) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
mod.tree_adult11 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + 
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult11)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult11: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree_adult11:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree_adult8:     A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult1:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.tree_adult11 12 5276.4 5346.8 -2626.2   5252.4              
## mod.tree_adult8  13 5269.0 5345.2 -2621.5   5243.0 9.4309      1
## mod.tree_adult1  14 5269.1 5351.2 -2620.6   5241.1 1.8897      1
##                  Pr(>Chisq)   
## mod.tree_adult11              
## mod.tree_adult8    0.002134 **
## mod.tree_adult1    0.169234   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.tree_adult12 <- glmer(cbind(survs, deaths) ~ site + pesticide*exclosure*census + 
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.tree_adult13 <- glmer(cbind(survs, deaths) ~ site + (pesticide+exclosure)*census + 
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.tree_adult14 <- glmer(cbind(survs, deaths) ~ site + pesticide*census + exclosure+ 
                           A.con_curt + A.het_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.tree_adult15 <- glmer(cbind(survs, deaths) ~ pesticide*census + exclosure+ 
                           A.con_curt + A.het_curt +
                     (1|plot/quadrat) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult13, mod.tree_adult14, mod.tree_adult15)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree_adult8:     A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree_adult1:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult15: cbind(survs, deaths) ~ pesticide * census + exclosure + A.con_curt + 
## mod.tree_adult15:     A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree_adult14: cbind(survs, deaths) ~ site + pesticide * census + exclosure + 
## mod.tree_adult14:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.tree_adult13: cbind(survs, deaths) ~ site + (pesticide + exclosure) * census + 
## mod.tree_adult13:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                  Df    AIC    BIC  logLik deviance   Chisq Chi Df
## mod.tree_adult8  13 5269.0 5345.2 -2621.5   5243.0               
## mod.tree_adult1  14 5269.1 5351.2 -2620.6   5241.1  1.8897      1
## mod.tree_adult15 22 5251.3 5380.2 -2603.7   5207.3 33.7936      8
## mod.tree_adult14 23 5251.5 5386.3 -2602.8   5205.5  1.8032      1
## mod.tree_adult13 26 5254.2 5406.6 -2601.1   5202.2  3.2734      3
##                  Pr(>Chisq)    
## mod.tree_adult8                
## mod.tree_adult1      0.1692    
## mod.tree_adult15  4.428e-05 ***
## mod.tree_adult14     0.1793    
## mod.tree_adult13     0.3514    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.tree_adult14, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide * census + exclosure +  
##     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5251.5   5386.3  -2602.8   5205.5     2568 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6563 -0.5954  0.3777  0.6731  2.3912 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2069   0.4548  
##  sp.         (Intercept) 0.5349   0.7314  
## Number of obs: 2591, groups:  quad.unique, 239; sp., 20
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -0.53085    0.70253  -0.756 0.449879    
## site1                  0.01310    0.06051   0.216 0.828655    
## site2                 -0.17776    0.05713  -3.111 0.001862 ** 
## pesticideF             0.73671    0.15200   4.847 1.25e-06 ***
## pesticideI             0.11356    0.15168   0.749 0.454041    
## pesticideC             0.03071    0.15330   0.200 0.841225    
## census16sp            -1.08980    0.22560  -4.831 1.36e-06 ***
## census17fa             0.34061    0.13116   2.597 0.009405 ** 
## census17sp            -0.46765    0.12700  -3.682 0.000231 ***
## exclosure1             0.09676    0.08302   1.166 0.243816    
## A.con_curt            -0.03045    0.01342  -2.269 0.023293 *  
## A.het_curt             0.13529    0.03718   3.639 0.000274 ***
## pesticideF:census16sp -0.89672    0.31745  -2.825 0.004731 ** 
## pesticideI:census16sp -0.51417    0.33445  -1.537 0.124198    
## pesticideC:census16sp -0.04847    0.33054  -0.147 0.883426    
## pesticideF:census17fa -0.92802    0.17629  -5.264 1.41e-07 ***
## pesticideI:census17fa -0.24958    0.18305  -1.363 0.172740    
## pesticideC:census17fa -0.34055    0.18703  -1.821 0.068632 .  
## pesticideF:census17sp -0.61059    0.17587  -3.472 0.000517 ***
## pesticideI:census17sp -0.20690    0.18029  -1.148 0.251123    
## pesticideC:census17sp -0.22216    0.18457  -1.204 0.228699    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 21 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## select the approprate model and test the sensitivity
mod.tree_adult.diags <- augment(mod.tree_adult1)
qplot(data=mod.tree_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.tree_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
qqPlot(resid(mod.tree_adult1))

qqPlot(resid(mod.tree_adult14))

library(sjPlot)
## Visit http://strengejacke.de/sjPlot for package-vignettes.
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
sjp.lmer(mod.tree_adult1, type='fe')

sjp.lmer(mod.tree_adult1, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.tree_adult1, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5269.1   5351.2  -2620.6   5241.1     2577 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6535 -0.5888  0.3808  0.6746  2.3480 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2088   0.4570  
##  sp.         (Intercept) 0.5085   0.7131  
## Number of obs: 2591, groups:  quad.unique, 239; sp., 20
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.31548    0.69975  -0.451 0.652105    
## site1        0.01390    0.06058   0.229 0.818543    
## site2       -0.17552    0.05721  -3.068 0.002155 ** 
## pesticideF   0.19337    0.11067   1.747 0.080599 .  
## pesticideI  -0.04799    0.11446  -0.419 0.675018    
## pesticideC  -0.13903    0.11502  -1.209 0.226767    
## exclosure1   0.09986    0.08308   1.202 0.229380    
## census16sp  -1.45024    0.12070 -12.015  < 2e-16 ***
## census17fa  -0.06942    0.06558  -1.058 0.289829    
## census17sp  -0.73059    0.06433 -11.356  < 2e-16 ***
## A.con_curt  -0.03002    0.01341  -2.239 0.025128 *  
## A.het_curt   0.13482    0.03723   3.621 0.000293 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## shrub
mod.shrub_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.sum_curt +
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
## model without sum area of adult neigboring trees
mod.shrub_adult2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.shrub_adult3 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census +                     (1|plot/quadrat) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.shrub_adult4 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + A.sum_curt +
                     (1|plot/quadrat) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
mod.shrub_adult5 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.shrub_adult6 <- glmer(cbind(survs, deaths) ~ site + pesticide*census + exclosure + 
                     (1|quad.unique) + (1|sp.),
                   data=subset(pestdat.ag_adult, growth.form=='shrub'), family=binomial, 
                   control=glmerControl(optimizer="optimx", 
                                        optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
anova(mod.shrub_adult1, mod.shrub_adult3, mod.shrub_adult4, mod.shrub_adult5, mod.shrub_adult6)
## Data: subset(pestdat.ag_adult, growth.form == "shrub")
## Models:
## mod.shrub_adult3: cbind(survs, deaths) ~ pesticide + exclosure + census + (1 | 
## mod.shrub_adult3:     plot/quadrat) + (1 | sp.)
## mod.shrub_adult4: cbind(survs, deaths) ~ pesticide + exclosure + census + A.sum_curt + 
## mod.shrub_adult4:     (1 | plot/quadrat) + (1 | sp.)
## mod.shrub_adult5: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.shrub_adult5:     (1 | quad.unique) + (1 | sp.)
## mod.shrub_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.shrub_adult1:     A.sum_curt + (1 | quad.unique) + (1 | sp.)
## mod.shrub_adult6: cbind(survs, deaths) ~ site + pesticide * census + exclosure + 
## mod.shrub_adult6:     (1 | quad.unique) + (1 | sp.)
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.shrub_adult3 11 486.90 546.98 -232.45   464.90              
## mod.shrub_adult4 12 488.18 553.73 -232.09   464.18 0.7110      1
## mod.shrub_adult5 12 486.44 551.98 -231.22   462.44 1.7487      0
## mod.shrub_adult1 13 487.26 558.27 -230.63   461.26 1.1793      1
## mod.shrub_adult6 21 496.75 611.46 -227.38   454.75 6.5059      8
##                  Pr(>Chisq)    
## mod.shrub_adult3               
## mod.shrub_adult4     0.3991    
## mod.shrub_adult5     <2e-16 ***
## mod.shrub_adult1     0.2775    
## mod.shrub_adult6     0.5908    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# add sum of neighoring trees make no sense

# model sensitivity tests
mod.shrub_adult.diags <- augment(mod.shrub_adult5)
qplot(data=mod.shrub_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.shrub_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.shrub_adult5))

sjp.lmer(mod.shrub_adult5, type='fe')

sjp.lmer(mod.shrub_adult5, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.shrub_adult5, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     (1 | quad.unique) + (1 | sp.)
##    Data: subset(pestdat.ag_adult, growth.form == "shrub")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    486.4    552.0   -231.2    462.4     1729 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -12.2932   0.0906   0.1282   0.1877   0.5913 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.35332  0.5944  
##  sp.         (Intercept) 0.09279  0.3046  
## Number of obs: 1741, groups:  quad.unique, 195; sp., 19
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.86943    0.56128   8.676  < 2e-16 ***
## site1       -0.15280    0.20412  -0.749 0.454108    
## site2        0.33838    0.23371   1.448 0.147657    
## pesticideF  -0.69626    0.40450  -1.721 0.085197 .  
## pesticideI  -0.17234    0.44019  -0.392 0.695412    
## pesticideC   0.02797    0.45794   0.061 0.951289    
## exclosure1   0.56311    0.31036   1.814 0.069616 .  
## census16sp  -1.65148    0.45466  -3.632 0.000281 ***
## census17fa   0.06338    0.58204   0.109 0.913291    
## census17sp  -1.11219    0.47702  -2.332 0.019724 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robi, I have one question here. For the tree survival analysis, the model considering interation between pesticide and census fits better than those without the interaction. Should we keep the interacton term? (See code: anova(mod.tree_adult1, mod.tree_adult8, mod.tree_adult13, mod.tree_adult14, mod.tree_adult15) Line 298) If so, how to interpret the results? I mean, several interactions between fungicide and censuses are negative. Is that mean they compared with the water treatment in first census?

Density dependence

##------------------------- Density dependence-------------------------##

## For all tree species together
ggplot(data=subset(pestdat.ag_adult, growth.form=='tree'), 
       aes(x=total, y=survs/total, colour=pesticide)) + 
  geom_jitter()+
  geom_smooth(aes(weight = total), method='glm', method.args=list(family="binomial")) +
  facet_wrap(~exclosure)

# hard to transform the density of seendlings to normal distributed
hist(pestdat.ag_adult$total)

hist(log10(pestdat.ag_adult$total))

pestdat.ag_adult$dens <- (pestdat.ag_adult$total - mean(pestdat.ag_adult$total))/(2*sd(pestdat.ag_adult$total))


mod.tree.ddst_adult1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt + dens + (1|quad.unique) + (1|sp.),
                             data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx",
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.tree.ddst_adult2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + A.het_curt + dens + (1|plot/quadrat) + (1|sp.),
                             data=subset(pestdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx",
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.tree_adult8, mod.tree.ddst_adult1, mod.tree.ddst_adult2)
## Data: subset(pestdat.ag_adult, growth.form == "tree")
## Models:
## mod.tree_adult8: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree_adult8:     A.het_curt + (1 | plot/quadrat) + (1 | sp.)
## mod.tree.ddst_adult2: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.tree.ddst_adult2:     A.het_curt + dens + (1 | plot/quadrat) + (1 | sp.)
## mod.tree.ddst_adult1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.tree.ddst_adult1:     A.con_curt + A.het_curt + dens + (1 | quad.unique) + (1 | 
## mod.tree.ddst_adult1:     sp.)
##                      Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.tree_adult8      13 5269.0 5345.2 -2621.5   5243.0              
## mod.tree.ddst_adult2 14 5265.7 5347.7 -2618.8   5237.7 5.3047      1
## mod.tree.ddst_adult1 15 5266.2 5354.1 -2618.1   5236.2 1.4530      1
##                      Pr(>Chisq)  
## mod.tree_adult8                  
## mod.tree.ddst_adult2    0.02127 *
## mod.tree.ddst_adult1    0.22805  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.tree.ddst_adult.diags <- augment(mod.tree.ddst_adult1)
qplot(data=mod.tree.ddst_adult.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data =mod.tree.ddst_adult.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.tree.ddst_adult1))

sjp.lmer(mod.tree.ddst_adult1, type='fe')

sjp.lmer(mod.tree.ddst_adult1, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.tree.ddst_adult1, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique) + (1 |      sp.)
##    Data: subset(pestdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5266.2   5354.1  -2618.1   5236.2     2576 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2386 -0.6019  0.3719  0.6590  2.3434 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2089   0.4570  
##  sp.         (Intercept) 0.4982   0.7058  
## Number of obs: 2591, groups:  quad.unique, 239; sp., 20
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.370323   0.700573  -0.529  0.59708    
## site1        0.009565   0.060692   0.158  0.87477    
## site2       -0.155523   0.057960  -2.683  0.00729 ** 
## pesticideF   0.212304   0.111082   1.911  0.05597 .  
## pesticideI  -0.049563   0.114523  -0.433  0.66517    
## pesticideC  -0.148329   0.115166  -1.288  0.19776    
## exclosure1   0.104011   0.083177   1.250  0.21113    
## census16sp  -1.516264   0.124482 -12.181  < 2e-16 ***
## census17fa  -0.083965   0.066092  -1.270  0.20393    
## census17sp  -0.757571   0.065556 -11.556  < 2e-16 ***
## A.con_curt  -0.027493   0.013494  -2.037  0.04161 *  
## A.het_curt   0.138247   0.037296   3.707  0.00021 ***
## dens        -0.091498   0.041339  -2.213  0.02687 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
##------------------------- individual species -------------------------------##
## And now pulling the species apart. 
## First need to filter down to species that can be analysed - need enough
## reps (lets say 10) and variation in abundance (at least > 5 
sp.counts <- summarise(group_by(pestdat.ag_adult, sp.), 
                       n.quadrat=length(total), max.dens=max(total), min.dens=min(total))

sp.counts$range <- sp.counts$max.dens - sp.counts$min.dens
sp.sel <- sp.counts$sp.[sp.counts$n.quadrat > 10 & sp.counts$range > 5]
sp.counts[sp.counts$sp. %in% sp.sel,]
## # A tibble: 8 × 5
##      sp. n.quadrat max.dens min.dens range
##   <fctr>     <int>    <int>    <int> <int>
## 1   ABNE        46       12        1    11
## 2   ACBA       260       19        1    18
## 3   ACPS       296       12        1    11
## 4   ACTE       147        8        1     7
## 5   FRMA       765       23        1    22
## 6   PIKO        75        8        1     7
## 7   QUMO        52        8        1     7
## 8   TIAM       673       22        1    21
ggplot(data=subset(pestdat.ag_adult, sp. %in% sp.sel), 
       aes(x=total, y=survs/total, colour=pesticide)) + 
  geom_jitter()+
  geom_smooth(aes(weight = total), method='glm',
              method.args=list(family="binomial"), fullrange=T) +
  facet_wrap(~ sp.,  scales='free')

# it seems only eight species have enough data, but conference intervals of ACMO and QUMO are too huge

## quick plot of only the species for which there appear to be enough data
# sp.sel2 <- c('ABNE', 'ACBA', 'ACMO', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'QUMO', 'TIAM')
sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')

ggplot(data=subset(pestdat.ag_adult, sp. %in% sp.sel2), 
       aes(x=total, y=survs/total, colour=pesticide)) + 
  geom_jitter()+
  geom_smooth(aes(weight = total), method='glm',
              method.args=list(family="binomial"), fullrange=T) +
  facet_wrap(~ sp. + exclosure,  scales='free')

nddmods <- sapply(sp.sel2, function(sp){
  print(sp)
  spdat <- filter(pestdat.ag_adult, sp. == sp)
  ## to increase ability to fit models, standardise total
  spdat$dens <- (spdat$total - mean(spdat$total))/(2*sd(spdat$total))
  spdat <- droplevels(spdat)
  mod <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census +
                 A.con_curt + A.het_curt + dens + (1|quad.unique), 
               data=spdat, family=binomial, 
               control=glmerControl(optimizer="optimx",
                                    optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
  )
}, simplify=FALSE)
## [1] "ACBA"
## [1] "ACPS"
## [1] "ACTE"
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] "FRMA"
## [1] "PIKO"
## [1] "TIAM"
## diagnostic plots
indmods.resids <- lapply(nddmods, augment) 

indmods.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.resids)), dat=indmods.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'

library(lattice)
lapply(nddmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACBA
## $ACBA$quad.unique

## 
## 
## $ACPS
## $ACPS$quad.unique

## 
## 
## $ACTE
## $ACTE$quad.unique

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## 
## $PIKO
## $PIKO$quad.unique

## 
## 
## $TIAM
## $TIAM$quad.unique

lapply(nddmods, function(x) summary(x, correlation=FALSE))
## $ACBA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    383.7    433.5   -177.8    355.7      246 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9506 -0.6146  0.3482  0.5208  1.3918 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 1.378    1.174   
## Number of obs: 260, groups:  quad.unique, 114
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  3.43543    4.08147   0.842  0.39995   
## site2       -1.36234    0.51243  -2.659  0.00785 **
## site3       -1.04899    0.61074  -1.718  0.08588 . 
## pesticideF  -0.27307    0.52806  -0.517  0.60508   
## pesticideI  -0.38875    0.52778  -0.737  0.46138   
## pesticideC  -0.05273    0.56505  -0.093  0.92565   
## exclosure1  -0.59540    0.43482  -1.369  0.17091   
## census16sp  -0.70149    0.60670  -1.156  0.24759   
## census17fa   0.42763    0.36927   1.158  0.24685   
## census17sp  -1.01130    0.35893  -2.817  0.00484 **
## A.con_curt  -0.01353    0.32369  -0.042  0.96667   
## A.het_curt  -0.04527    0.22333  -0.203  0.83937   
## dens        -0.44753    0.28677  -1.561  0.11861   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE)  or
##   vcov(....)  if you need it
## 
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    370.1    421.7   -171.0    342.1      282 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.3717  0.0843  0.2185  0.4799  1.1045 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.4266   0.6531  
## Number of obs: 296, groups:  quad.unique, 159
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  4.60321    2.84794   1.616  0.10602   
## site2       -0.08561    0.31637  -0.271  0.78671   
## site3        0.31949    0.41168   0.776  0.43772   
## pesticideF   0.08305    0.35728   0.232  0.81618   
## pesticideI   0.12886    0.38874   0.332  0.74029   
## pesticideC  -0.08517    0.39967  -0.213  0.83125   
## exclosure1  -0.22560    0.33685  -0.670  0.50302   
## census16sp  -2.49543    1.27919  -1.951  0.05108 . 
## census17fa  -3.20721    1.02650  -3.124  0.00178 **
## census17sp  -1.67945    1.10957  -1.514  0.13012   
## A.con_curt  -0.11535    0.15830  -0.729  0.46617   
## A.het_curt   0.04720    0.15061   0.313  0.75397   
## dens         0.22956    0.22045   1.041  0.29772   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE)  or
##   vcov(....)  if you need it
## 
## $ACTE
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    188.1    227.0    -81.1    162.1      134 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8905 -0.1223  0.3494  0.5458  1.2432 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0        0       
## Number of obs: 147, groups:  quad.unique, 59
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -5.62217    4.29189  -1.310   0.1902  
## site2       -0.15382    0.43979  -0.350   0.7265  
## pesticideF  -0.73375    0.55433  -1.324   0.1856  
## pesticideI  -0.63145    0.52160  -1.211   0.2260  
## pesticideC  -0.51642    0.56344  -0.917   0.3594  
## exclosure1  -0.11001    0.55640  -0.198   0.8433  
## census16sp  -1.48944    0.62940  -2.366   0.0180 *
## census17fa   0.44397    0.45710   0.971   0.3314  
## census17sp   0.20571    0.47039   0.437   0.6619  
## A.con_curt  -0.47809    0.28447  -1.681   0.0928 .
## A.het_curt   0.54951    0.23876   2.301   0.0214 *
## dens        -0.05784    0.27089  -0.214   0.8309  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1680.9   1745.8   -826.4   1652.9      751 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9632 -0.6714  0.3874  0.6946  1.9496 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2132   0.4617  
## Number of obs: 765, groups:  quad.unique, 230
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.050470   0.929449   0.054   0.9567    
## site2       -0.354717   0.175897  -2.017   0.0437 *  
## site3       -0.180762   0.179258  -1.008   0.3133    
## pesticideF   0.062224   0.163599   0.380   0.7037    
## pesticideI  -0.301626   0.162044  -1.861   0.0627 .  
## pesticideC  -0.139590   0.165317  -0.844   0.3985    
## exclosure1   0.251150   0.120128   2.091   0.0366 *  
## census16sp  -1.560813   0.176617  -8.837  < 2e-16 ***
## census17fa   0.310369   0.136724   2.270   0.0232 *  
## census17sp  -0.495096   0.108348  -4.570 4.89e-06 ***
## A.con_curt   0.006949   0.016989   0.409   0.6825    
## A.het_curt   0.114272   0.052438   2.179   0.0293 *  
## dens        -0.179276   0.120689  -1.485   0.1374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE)  or
##   vcov(....)  if you need it
## 
## $PIKO
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    146.7    179.2    -59.4    118.7       61 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8041 -0.5187  0.3436  0.6499  1.1901 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.6002   0.7748  
## Number of obs: 75, groups:  quad.unique, 43
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -6.89039    6.98760  -0.986  0.32409   
## site2       -0.17103    0.83838  -0.204  0.83836   
## site3       -0.80554    0.95573  -0.843  0.39931   
## pesticideF   2.61863    0.95197   2.751  0.00595 **
## pesticideI   0.66617    0.85409   0.780  0.43541   
## pesticideC   0.76070    0.99930   0.761  0.44651   
## exclosure1   0.24621    0.62163   0.396  0.69206   
## census16sp  -1.45947    1.96659  -0.742  0.45801   
## census17fa   0.53421    0.71735   0.745  0.45645   
## census17sp  -0.46806    0.59740  -0.784  0.43334   
## A.con_curt   0.62626    0.42384   1.478  0.13952   
## A.het_curt   0.05408    0.24498   0.221  0.82529   
## dens         0.07058    0.51695   0.136  0.89140   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE)  or
##   vcov(....)  if you need it
## 
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + dens + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1992.2   2055.4   -982.1   1964.2      659 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4000 -0.6972  0.0891  0.8041  2.1376 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2655   0.5153  
## Number of obs: 673, groups:  quad.unique, 233
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.66161    1.02028   1.629  0.10340    
## site2       -0.16619    0.14091  -1.179  0.23826    
## site3        0.23276    0.14022   1.660  0.09693 .  
## pesticideF   0.38856    0.14601   2.661  0.00779 ** 
## pesticideI   0.14049    0.15422   0.911  0.36231    
## pesticideC  -0.18497    0.15618  -1.184  0.23628    
## exclosure1   0.16308    0.11200   1.456  0.14536    
## census16sp  -1.40203    0.46003  -3.048  0.00231 ** 
## census17fa  -0.25289    0.09225  -2.741  0.00612 ** 
## census17sp  -1.17234    0.10043 -11.673  < 2e-16 ***
## A.con_curt  -0.20009    0.04056  -4.933  8.1e-07 ***
## A.het_curt   0.08584    0.05139   1.670  0.09483 .  
## dens        -0.26475    0.08656  -3.059  0.00222 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(...., correlation=TRUE)  or
##   vcov(....)  if you need it
lapply(nddmods, function(x) anova(x))
## $ACBA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  2.6593  1.3296  1.3296
## pesticide   3  0.6135  0.2045  0.2045
## exclosure   1  1.0220  1.0220  1.0220
## census      3 14.4767  4.8256  4.8256
## A.con_curt  1  0.0015  0.0015  0.0015
## A.het_curt  1  0.0424  0.0424  0.0424
## dens        1  2.9350  2.9350  2.9350
## 
## $ACPS
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  0.6643  0.3322  0.3322
## pesticide   3  0.3228  0.1076  0.1076
## exclosure   1  0.8021  0.8021  0.8021
## census      3 15.2684  5.0895  5.0895
## A.con_curt  1  0.4563  0.4563  0.4563
## A.het_curt  1  0.0728  0.0728  0.0728
## dens        1  1.0919  1.0919  1.0919
## 
## $ACTE
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        1 0.2896  0.2896  0.2896
## pesticide   3 2.5762  0.8587  0.8587
## exclosure   1 0.1924  0.1924  0.1924
## census      3 8.5964  2.8655  2.8655
## A.con_curt  1 3.5758  3.5758  3.5758
## A.het_curt  1 5.2620  5.2620  5.2620
## dens        1 0.0456  0.0456  0.0456
## 
## $FRMA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  10.127   5.063  5.0634
## pesticide   3   4.128   1.376  1.3759
## exclosure   1   4.588   4.588  4.5884
## census      3 139.199  46.400 46.3997
## A.con_curt  1   0.326   0.326  0.3261
## A.het_curt  1   3.534   3.534  3.5341
## dens        1   2.478   2.478  2.4781
## 
## $PIKO
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 0.0797  0.0399  0.0399
## pesticide   3 9.7310  3.2437  3.2437
## exclosure   1 0.0473  0.0473  0.0473
## census      3 3.1894  1.0631  1.0631
## A.con_curt  1 2.1612  2.1612  2.1612
## A.het_curt  1 0.0727  0.0727  0.0727
## dens        1 0.0192  0.0192  0.0192
## 
## $TIAM
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  10.486   5.243  5.2429
## pesticide   3   9.284   3.095  3.0947
## exclosure   1   3.021   3.021  3.0215
## census      3 143.101  47.700 47.7003
## A.con_curt  1  33.836  33.836 33.8358
## A.het_curt  1   2.565   2.565  2.5647
## dens        1   9.901   9.901  9.9010
##---------------- Survival of other tree seedlings-------------------##
# sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')
## other tree species
pestdat.tree_adult1 <- pestdat.ag_adult[pestdat.ag_adult$growth.form == 'tree',][!(pestdat.ag_adult[pestdat.ag_adult$growth.form == 'tree',]$sp. %in% sp.sel2),]
  
summary(pestdat.tree_adult1)
##  site       quadrat           sp.      pesticide exclosure growth.form
##  1:139   Min.   :101.0   ACMO   :139   W: 83     0:172     shrub:  0  
##  2:198   1st Qu.:206.0   ACMA   : 80   F:126     1:203     tree :375  
##  3: 38   Median :305.0   QUMO   : 52   I: 97                          
##          Mean   :315.8   ABNE   : 46   C: 69                          
##          3rd Qu.:410.0   MAAM   : 12                                  
##          Max.   :510.0   ACTR   : 11                                  
##                          (Other): 35                                  
##   census        deaths           total            survs      
##  16fa: 76   Min.   :0.0000   Min.   : 1.000   Min.   :0.000  
##  16sp: 37   1st Qu.:0.0000   1st Qu.: 1.000   1st Qu.:1.000  
##  17fa:189   Median :0.0000   Median : 1.000   Median :1.000  
##  17sp: 73   Mean   :0.3387   Mean   : 1.555   Mean   :1.216  
##             3rd Qu.:0.0000   3rd Qu.: 1.000   3rd Qu.:1.000  
##             Max.   :6.0000   Max.   :12.000   Max.   :9.000  
##                                                              
##   quad.unique   plot             quad.sp        A.sum     
##  2_1_209: 15   1_0: 50   1_0_105_ACTR:  4   Min.   :4071  
##  2_0_206:  8   1_1: 89   1_0_406_ACMA:  4   1st Qu.:4906  
##  1_1_407:  7   2_0:103   1_1_305_ACMA:  4   Median :5403  
##  1_1_410:  7   2_1: 95   1_1_308_ACMO:  4   Mean   :5562  
##  2_0_204:  7   3_0: 19   1_1_406_ACMO:  4   3rd Qu.:6041  
##  2_1_206:  7   3_1: 19   1_1_407_ACMO:  4   Max.   :8723  
##  (Other):324             (Other)     :351                 
##      A.con              A.het         con.dens        A.con_curt    
##  Min.   :   0.000   Min.   :2865   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:   8.913   1st Qu.:4499   1st Qu.: 1.000   1st Qu.: 2.073  
##  Median : 206.523   Median :5105   Median : 6.000   Median : 5.911  
##  Mean   : 344.675   Mean   :5218   Mean   : 7.691   Mean   : 5.172  
##  3rd Qu.: 358.975   3rd Qu.:5694   3rd Qu.: 9.000   3rd Qu.: 7.107  
##  Max.   :3114.708   Max.   :8560   Max.   :46.000   Max.   :14.604  
##                                                                     
##    A.het_curt      A.sum_curt    con.dens_curt        dens        
##  Min.   :14.20   Min.   :15.97   Min.   :0.000   Min.   :-0.3037  
##  1st Qu.:16.51   1st Qu.:16.99   1st Qu.:1.000   1st Qu.:-0.3037  
##  Median :17.22   Median :17.55   Median :1.817   Median :-0.3037  
##  Mean   :17.27   Mean   :17.67   Mean   :1.510   Mean   :-0.2098  
##  3rd Qu.:17.86   3rd Qu.:18.21   3rd Qu.:2.080   3rd Qu.:-0.3037  
##  Max.   :20.46   Max.   :20.59   Max.   :3.583   Max.   : 1.5574  
## 
mod.s.t_adult.r1 <- glmer(cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
                                 A.con_curt + A.het_curt + (1|quad.unique) + (1|sp.),
                             data=pestdat.tree_adult1, family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.s.t_adult.r2 <- glmer(cbind(survs, deaths) ~ pesticide + exclosure + census + 
                                 A.con_curt + A.het_curt + (1|plot/quad.unique) + (1|sp.),
                             data=pestdat.tree_adult1, family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.046875 (tol =
## 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
anova(mod.s.t_adult.r1, mod.s.t_adult.r2)
## Data: pestdat.tree_adult1
## Models:
## mod.s.t_adult.r2: cbind(survs, deaths) ~ pesticide + exclosure + census + A.con_curt + 
## mod.s.t_adult.r2:     A.het_curt + (1 | plot/quad.unique) + (1 | sp.)
## mod.s.t_adult.r1: cbind(survs, deaths) ~ site + pesticide + exclosure + census + 
## mod.s.t_adult.r1:     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.s.t_adult.r2 13 421.61 472.66 -197.81   395.61              
## mod.s.t_adult.r1 14 422.78 477.76 -197.39   394.78 0.8275      1
##                  Pr(>Chisq)
## mod.s.t_adult.r2           
## mod.s.t_adult.r1      0.363
# first, check the random effects
sjp.lmer(mod.s.t_adult.r2, type='re', sort.est="sort.all", p.kr = FALSE)
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## setting plot as random effects doesn't take any variance, the first model with site as fixed effect woulb be appropriate

mod.s.t_adult.r.diags <- augment(mod.s.t_adult.r1)
qplot(data=mod.s.t_adult.r.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data=mod.s.t_adult.r.diags, x=pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

## plot fixed and random effects
sjp.lmer(mod.s.t_adult.r1, type='fe', p.kr = FALSE)

sjp.lmer(mod.s.t_adult.r1, type='re', sort.est="sort.all", p.kr = FALSE)
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.s.t_adult.r1))

summary(mod.s.t_adult.r1, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + pesticide + exclosure + census +  
##     A.con_curt + A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: pestdat.tree_adult1
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    422.8    477.8   -197.4    394.8      361 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7012  0.2109  0.3113  0.4231  1.4231 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.4184   0.6469  
##  sp.         (Intercept) 0.9803   0.9901  
## Number of obs: 375, groups:  quad.unique, 138; sp., 14
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  4.38646    2.71314   1.617  0.10593   
## site1       -0.05411    0.27860  -0.194  0.84601   
## site2       -0.23912    0.28191  -0.848  0.39632   
## pesticideF   0.27889    0.40408   0.690  0.49007   
## pesticideI   0.09446    0.39454   0.239  0.81079   
## pesticideC  -0.24100    0.41618  -0.579  0.56253   
## exclosure1  -0.31163    0.30125  -1.034  0.30091   
## census16sp  -1.64748    0.56673  -2.907  0.00365 **
## census17fa  -0.28476    0.50103  -0.568  0.56981   
## census17sp  -0.97794    0.51302  -1.906  0.05662 . 
## A.con_curt  -0.07461    0.07035  -1.060  0.28890   
## A.het_curt  -0.08986    0.14921  -0.602  0.54702   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# no significant effect for the rest tree seedlings

Overall, fungicide slightly higher survival of tree seedlings pooled. Effect of conspecific neighboring trees is negative while effect of heterospecific neighboring trees is positive. Seedling density shows negative density dependence. When tested at species level. Results are consistent with that at community-level.Fungicide heigher survival of two out of nine abundant species (PIKO and TIAM). Conspecific neighboring trees have negative effect on seedling survival of TIAM. Heterospecific neighboring trees have positive effect on ACTE and FRAM. Negative density shows in FRMA and TIAM. Fungicide heigher suvival when the rest of tree species pooled (i.e., not including the common species).

Analysis of Tree Recruitment

Tree recruitments mainly occur in the first census (spring), we only use two spring censuses (16sp and 17sp) in testing recruitment.

###-------------------- recruitment annlysis -----------------------------###
pest.rect.ag_adult <- read.csv("data/pest.rect.ag_adult.csv")
pest.rect.ag_adult$pesticide <- factor(pest.rect.ag_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pest.rect.ag_adult$census <- factor(pest.rect.ag_adult$census)
pest.rect.ag_adult$exclosure <- factor(pest.rect.ag_adult$exclosure)
pest.rect.ag_adult$site <- factor(pest.rect.ag_adult$site)
summary(pest.rect.ag_adult)
##        X        site       quadrat           sp.      pesticide exclosure
##  Min.   :   1   1:497   Min.   :101.0   TIAM   :410   W:352     0:715    
##  1st Qu.: 357   2:564   1st Qu.:204.0   FRMA   :283   F:401     1:710    
##  Median : 713   3:364   Median :304.0   ACPS   :198   I:357              
##  Mean   : 713           Mean   :304.9   ACBA   :126   C:315              
##  3rd Qu.:1069           3rd Qu.:408.0   ACMO   : 70                      
##  Max.   :1425           Max.   :510.0   ACTE   : 64                      
##                                         (Other):274                      
##  growth.form   census       recruits       quad.unique    plot    
##  shrub: 113   16sp:747   Min.   : 1.000   2_1_302:  13   1_0:228  
##  tree :1312   17sp:678   1st Qu.: 1.000   2_0_504:  12   1_1:269  
##                          Median : 2.000   2_0_505:  12   2_0:303  
##                          Mean   : 3.331   1_1_410:  11   2_1:261  
##                          3rd Qu.: 5.000   2_0_310:  11   3_0:184  
##                          Max.   :21.000   1_0_410:  10   3_1:180  
##                                           (Other):1356            
##          quad.sp         A.sum           A.con             A.het      
##  1_0_101_ACPS:   2   Min.   : 3990   Min.   :   0.00   Min.   : 2374  
##  1_0_101_TIAM:   2   1st Qu.: 4991   1st Qu.:  75.49   1st Qu.: 4039  
##  1_0_102_FRMA:   2   Median : 5463   Median : 398.78   Median : 4877  
##  1_0_102_TIAM:   2   Mean   : 5639   Mean   : 742.11   Mean   : 4897  
##  1_0_103_ACMO:   2   3rd Qu.: 6084   3rd Qu.:1248.58   3rd Qu.: 5573  
##  1_0_103_TIAM:   2   Max.   :10411   Max.   :5308.94   Max.   :10411  
##  (Other)     :1413                                                    
##     con.dens       A.con_curt       A.het_curt      A.sum_curt   
##  Min.   : 0.00   Min.   : 0.000   Min.   :13.34   Min.   :15.86  
##  1st Qu.: 1.00   1st Qu.: 4.226   1st Qu.:15.93   1st Qu.:17.09  
##  Median : 8.00   Median : 7.361   Median :16.96   Median :17.61  
##  Mean   :10.22   Mean   : 7.087   Mean   :16.87   Mean   :17.74  
##  3rd Qu.:14.00   3rd Qu.:10.768   3rd Qu.:17.73   3rd Qu.:18.26  
##  Max.   :46.00   Max.   :17.445   Max.   :21.84   Max.   :21.84  
##                                                                  
##  con.dens_curt  
##  Min.   :0.000  
##  1st Qu.:1.000  
##  Median :2.000  
##  Mean   :1.766  
##  3rd Qu.:2.410  
##  Max.   :3.583  
## 
###------------- Recruitment models ------------------###
pest.rect.ag_adult$site <- factor(pest.rect.ag_adult$site) 
contrasts(pest.rect.ag_adult$site) <- "contr.sum"

## Only tree species
pest.t.rect_adult <- pest.rect.ag_adult[pest.rect.ag_adult$growth.form == 'tree',]
mod.pest.recr1 <- glmer(recruits ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
                          (1|quad.unique) + (1|sp.),
                        data=pest.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

sjp.lmer(mod.pest.recr1, type='fe')

sjp.lmer(mod.pest.recr1, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Random effects of FRMA and TIAM go extramely, try to put group them against other species as the fixed effect
pest.t.rect_adult$sp.group <- ifelse(pest.t.rect_adult$sp.=='FRMA' | pest.t.rect_adult$sp.=='TIAM', 1, 0)

pest.t.rect_adult$sp.group <- factor(pest.t.rect_adult$sp.group) 
contrasts(pest.t.rect_adult$site) <- "contr.sum"
summary(pest.t.rect_adult)
##        X          site       quadrat         sp.      pesticide exclosure
##  Min.   :   1.0   1:452   Min.   :101   TIAM   :410   W:325     0:657    
##  1st Qu.: 330.8   2:516   1st Qu.:204   FRMA   :283   F:368     1:655    
##  Median : 658.5   3:344   Median :304   ACPS   :198   I:327              
##  Mean   : 664.2           Mean   :305   ACBA   :126   C:292              
##  3rd Qu.: 986.2           3rd Qu.:408   ACMO   : 70                      
##  Max.   :1425.0           Max.   :510   ACTE   : 64                      
##                                         (Other):161                      
##  growth.form   census       recruits       quad.unique    plot    
##  shrub:   0   16sp:689   Min.   : 1.000   2_1_302:  12   1_0:208  
##  tree :1312   17sp:623   1st Qu.: 1.000   2_0_504:  11   1_1:244  
##                          Median : 2.000   1_1_207:  10   2_0:276  
##                          Mean   : 3.528   1_1_503:  10   2_1:240  
##                          3rd Qu.: 5.000   2_0_505:  10   3_0:173  
##                          Max.   :21.000   2_1_209:  10   3_1:171  
##                                           (Other):1249            
##          quad.sp         A.sum           A.con            A.het      
##  1_0_101_ACPS:   2   Min.   : 3990   Min.   :   0.0   Min.   : 2374  
##  1_0_101_TIAM:   2   1st Qu.: 4991   1st Qu.: 159.0   1st Qu.: 3942  
##  1_0_102_FRMA:   2   Median : 5463   Median : 507.2   Median : 4805  
##  1_0_102_TIAM:   2   Mean   : 5636   Mean   : 806.0   Mean   : 4830  
##  1_0_103_ACMO:   2   3rd Qu.: 6086   3rd Qu.:1345.7   3rd Qu.: 5496  
##  1_0_103_TIAM:   2   Max.   :10411   Max.   :5308.9   Max.   :10186  
##  (Other)     :1300                                                   
##     con.dens      A.con_curt       A.het_curt      A.sum_curt   
##  Min.   : 0.0   Min.   : 0.000   Min.   :13.34   Min.   :15.86  
##  1st Qu.: 3.0   1st Qu.: 5.417   1st Qu.:15.80   1st Qu.:17.09  
##  Median : 9.0   Median : 7.975   Median :16.87   Median :17.61  
##  Mean   :11.1   Mean   : 7.692   Mean   :16.79   Mean   :17.74  
##  3rd Qu.:14.0   3rd Qu.:11.040   3rd Qu.:17.65   3rd Qu.:18.26  
##  Max.   :46.0   Max.   :17.445   Max.   :21.68   Max.   :21.84  
##                                                                 
##  con.dens_curt   sp.group
##  Min.   :0.000   0:619   
##  1st Qu.:1.442   1:693   
##  Median :2.080           
##  Mean   :1.915           
##  3rd Qu.:2.410           
##  Max.   :3.583           
## 
mod.pest.recr2 <- glmer(recruits ~ site + pesticide + exclosure + census + 
                          A.con_curt + A.het_curt + sp.group +
                          (1|quad.unique) + (1|sp.),
                        data=pest.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
mod.pest.recr3 <- glmer(recruits ~ pesticide + exclosure + census + 
                          A.con_curt + A.het_curt + sp.group +
                          (1|plot/quadrat) + (1|sp.),
                        data=pest.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.pest.recr1, mod.pest.recr2, mod.pest.recr3)
## Data: pest.t.rect_adult
## Models:
## mod.pest.recr1: recruits ~ site + pesticide + exclosure + census + A.con_curt + 
## mod.pest.recr1:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.recr3: recruits ~ pesticide + exclosure + census + A.con_curt + A.het_curt + 
## mod.pest.recr3:     sp.group + (1 | plot/quadrat) + (1 | sp.)
## mod.pest.recr2: recruits ~ site + pesticide + exclosure + census + A.con_curt + 
## mod.pest.recr2:     A.het_curt + sp.group + (1 | quad.unique) + (1 | sp.)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.pest.recr1 12 5921.6 5983.8 -2948.8   5897.6                         
## mod.pest.recr3 12 5918.0 5980.1 -2947.0   5894.0 3.6973      0    < 2e-16
## mod.pest.recr2 13 5914.3 5981.6 -2944.1   5888.3 5.6902      1    0.01706
##                   
## mod.pest.recr1    
## mod.pest.recr3 ***
## mod.pest.recr2 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.pest.recr.diags2 <- augment(mod.pest.recr2)
qplot(data=mod.pest.recr.diags2, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.pest.recr.diags2, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.pest.recr2))

sjp.lmer(mod.pest.recr2, type='fe')

sjp.lmer(mod.pest.recr2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.pest.recr2, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + sp.group + (1 | quad.unique) + (1 | sp.)
##    Data: pest.t.rect_adult
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5914.3   5981.6  -2944.1   5888.3     1299 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6577 -0.8584 -0.3143  0.4867  7.4046 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.06006  0.2451  
##  sp.         (Intercept) 0.08898  0.2983  
## Number of obs: 1312, groups:  quad.unique, 238; sp., 17
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.647454   0.385920   1.678 0.093407 .  
## site1       -0.038577   0.032841  -1.175 0.240136    
## site2        0.212409   0.031561   6.730  1.7e-11 ***
## pesticideF   0.204853   0.060803   3.369 0.000754 ***
## pesticideI   0.001089   0.063150   0.017 0.986246    
## pesticideC  -0.088692   0.064183  -1.382 0.167014    
## exclosure1   0.002293   0.045819   0.050 0.960095    
## census17sp  -0.397786   0.033850 -11.751  < 2e-16 ***
## A.con_curt   0.027352   0.007757   3.526 0.000422 ***
## A.het_curt  -0.011216   0.020534  -0.546 0.584906    
## sp.group1    0.889435   0.238743   3.725 0.000195 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
##The latter model considering fits better, but results do not's change, so it's ok to use the previous model 


###---------------- individual species analysis ----------------------------###
## reps (lets say 10) and variation in abundance (at least >5
sp.rect.counts <- summarise(group_by(pest.t.rect_adult, sp.), 
                       n.quadrat=length(recruits), max.rect=max(recruits), min.rect=min(recruits))

sp.rect.counts$range <- sp.rect.counts$max.rect - sp.rect.counts$min.rect
sp.rect.sel <- sp.rect.counts$sp.[sp.rect.counts$n.quadrat > 10 & sp.rect.counts$range > 5]
sp.rect.counts[sp.rect.counts$sp. %in% sp.rect.sel,]
## # A tibble: 8 × 5
##      sp. n.quadrat max.rect min.rect range
##   <fctr>     <int>    <int>    <int> <int>
## 1   ABNE        45       12        1    11
## 2   ACBA       126       13        1    12
## 3   ACPS       198       12        1    11
## 4   ACTE        64        7        1     6
## 5   FRMA       283       21        1    20
## 6   PIKO        43        8        1     7
## 7   QUMO        40        8        1     7
## 8   TIAM       410       21        1    20
# sp.rect.sel2 <- c('ABNE', 'ACBA',  'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')
sp.rect.sel2 <- c('ACBA',  'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')

rectmods <- sapply(sp.rect.sel2, function(sp){
  print(sp)
  spdat <- filter(pest.t.rect_adult, sp. == sp)
  spdat <- droplevels(spdat)
  mod <- glmer(recruits ~ site + pesticide + exclosure + census +
                 A.con_curt + A.het_curt + (1|quad.unique), 
               data=spdat, family=poisson, 
               control=glmerControl(optimizer="optimx",
                                    optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
  )
}, simplify=FALSE)
## [1] "ACBA"
## [1] "ACPS"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0323344 (tol =
## 0.001, component 1)
## [1] "ACTE"
## [1] "FRMA"
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## [1] "PIKO"
## [1] "QUMO"
## [1] "TIAM"
## diagnostic plots
indmods.rect.resids <- lapply(rectmods, augment) 

indmods.rect.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.rect.resids)), dat=indmods.rect.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.rect.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'

lapply(rectmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACBA
## $ACBA$quad.unique

## 
## 
## $ACPS
## $ACPS$quad.unique

## 
## 
## $ACTE
## $ACTE$quad.unique

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## 
## $PIKO
## $PIKO$quad.unique

## 
## 
## $QUMO
## $QUMO$quad.unique

## 
## 
## $TIAM
## $TIAM$quad.unique

lapply(rectmods, function(x) summary(x, correlation=FALSE))
## $ACBA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    406.5    437.7   -192.2    384.5      115 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8758 -0.4577 -0.2019  0.2797  2.7061 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1026   0.3203  
## Number of obs: 126, groups:  quad.unique, 101
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  1.65540    1.61951   1.022   0.3067  
## site2       -0.30498    0.19797  -1.540   0.1234  
## site3       -0.30455    0.24054  -1.266   0.2055  
## pesticideF  -0.13314    0.21391  -0.622   0.5337  
## pesticideI  -0.08852    0.20630  -0.429   0.6679  
## pesticideC  -0.42020    0.22896  -1.835   0.0665 .
## exclosure1  -0.33281    0.17976  -1.851   0.0641 .
## census17sp   0.02200    0.14881   0.148   0.8825  
## A.con_curt   0.06581    0.13318   0.494   0.6212  
## A.het_curt  -0.04269    0.08931  -0.478   0.6326  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    710.8    746.9   -344.4    688.8      187 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3584 -0.5834 -0.1737  0.3254  2.7500 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.08081  0.2843  
## Number of obs: 198, groups:  quad.unique, 157
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.547614   0.002990  183.16   <2e-16 ***
## site2        0.001170   0.002985    0.39    0.695    
## site3       -0.200816   0.002985  -67.27   <2e-16 ***
## pesticideF   0.367167   0.002986  122.97   <2e-16 ***
## pesticideI   0.200777   0.002985   67.25   <2e-16 ***
## pesticideC  -0.111387   0.002984  -37.32   <2e-16 ***
## exclosure1   0.406537   0.002988  136.05   <2e-16 ***
## census17sp   0.631224   0.002989  211.20   <2e-16 ***
## A.con_curt   0.028199   0.002926    9.64   <2e-16 ***
## A.het_curt  -0.043298   0.002237  -19.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0323344 (tol = 0.001, component 1)
## 
## 
## $ACTE
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    195.3    216.9    -87.7    175.3       54 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9890 -0.3944 -0.1977  0.1924  3.7139 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0        0       
## Number of obs: 64, groups:  quad.unique, 53
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.86641    2.09376   0.891    0.373
## site2        0.13401    0.22485   0.596    0.551
## pesticideF  -0.39642    0.30706  -1.291    0.197
## pesticideI  -0.17990    0.27248  -0.660    0.509
## pesticideC   0.00579    0.31724   0.018    0.985
## exclosure1  -0.19940    0.33783  -0.590    0.555
## census17sp   0.16891    0.21189   0.797    0.425
## A.con_curt   0.24979    0.16561   1.508    0.131
## A.het_curt  -0.13422    0.11351  -1.183    0.237
## 
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1274.0   1314.1   -626.0   1252.0      272 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.42007 -0.51942 -0.03114  0.42931  3.05340 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1421   0.377   
## Number of obs: 283, groups:  quad.unique, 225
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.1996672  0.6327376  -0.316   0.7523    
## site2        0.8794064  0.1028044   8.554  < 2e-16 ***
## site3        0.4617114  0.1104352   4.181  2.9e-05 ***
## pesticideF  -0.0304349  0.1090314  -0.279   0.7801    
## pesticideI   0.0100566  0.1095858   0.092   0.9269    
## pesticideC  -0.0768766  0.1105538  -0.695   0.4868    
## exclosure1   0.1397609  0.0801840   1.743   0.0813 .  
## census17sp  -1.5931387  0.1184738 -13.447  < 2e-16 ***
## A.con_curt   0.0003858  0.0115886   0.033   0.9734    
## A.het_curt   0.0738762  0.0357336   2.067   0.0387 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## 
## 
## $PIKO
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    168.8    188.1    -73.4    146.8       32 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0326 -0.4647 -0.2246  0.1686  2.0742 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1135   0.3369  
## Number of obs: 43, groups:  quad.unique, 42
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  4.55983    2.90398   1.570   0.1164  
## site2       -0.31798    0.33577  -0.947   0.3436  
## site3        0.24589    0.35496   0.693   0.4885  
## pesticideF  -0.02219    0.34243  -0.065   0.9483  
## pesticideI  -0.40482    0.37751  -1.072   0.2836  
## pesticideC  -0.24694    0.41616  -0.593   0.5529  
## exclosure1  -0.36151    0.26560  -1.361   0.1735  
## census17sp  -0.58356    0.28115  -2.076   0.0379 *
## A.con_curt  -0.25522    0.18845  -1.354   0.1756  
## A.het_curt  -0.04435    0.09835  -0.451   0.6520  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $QUMO
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    143.6    162.2    -60.8    121.6       29 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1079 -0.4019 -0.1627  0.2458  2.9716 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0        0       
## Number of obs: 40, groups:  quad.unique, 40
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -3.83219    4.02768  -0.952   0.3414  
## site2       -0.40959    0.62427  -0.656   0.5118  
## site3       -0.66936    0.36917  -1.813   0.0698 .
## pesticideF   0.20670    0.31899   0.648   0.5170  
## pesticideI   0.44270    0.30623   1.446   0.1483  
## pesticideC  -0.06282    0.42210  -0.149   0.8817  
## exclosure1  -0.33106    0.44759  -0.740   0.4595  
## census17sp  -0.16427    0.50083  -0.328   0.7429  
## A.con_curt   0.14229    0.11441   1.244   0.2136  
## A.het_curt   0.21224    0.21168   1.003   0.3160  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   2086.3   2130.5  -1032.1   2064.3      399 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4275 -0.8542 -0.2085  0.5745  5.1464 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1257   0.3545  
## Number of obs: 410, groups:  quad.unique, 232
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.87154    0.63050   1.382  0.16688    
## site2        0.07950    0.08648   0.919  0.35796    
## site3       -0.25182    0.08783  -2.867  0.00414 ** 
## pesticideF   0.38493    0.08891   4.329  1.5e-05 ***
## pesticideI  -0.01172    0.09504  -0.123  0.90187    
## pesticideC  -0.09726    0.09680  -1.005  0.31503    
## exclosure1  -0.18156    0.06896  -2.633  0.00847 ** 
## census17sp  -0.43175    0.04661  -9.263  < 2e-16 ***
## A.con_curt   0.04445    0.02466   1.802  0.07147 .  
## A.het_curt   0.02319    0.03190   0.727  0.46737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lapply(rectmods, function(x) anova(x))
## $ACBA
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 3.1207  1.5604  1.5604
## pesticide   3 3.0695  1.0232  1.0232
## exclosure   1 4.5944  4.5944  4.5944
## census      1 0.0575  0.0575  0.0575
## A.con_curt  1 0.3377  0.3377  0.3377
## A.het_curt  1 0.2235  0.2235  0.2235
## 
## $ACPS
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  2.8820  1.4410  1.4410
## pesticide   3 14.9578  4.9859  4.9859
## exclosure   1 16.0084 16.0084 16.0084
## census      1 23.6082 23.6082 23.6082
## A.con_curt  1  0.2093  0.2093  0.2093
## A.het_curt  1  0.5247  0.5247  0.5247
## 
## $ACTE
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        1 0.04332 0.04332  0.0433
## pesticide   3 2.14198 0.71399  0.7140
## exclosure   1 2.35134 2.35134  2.3513
## census      1 0.46046 0.46046  0.4605
## A.con_curt  1 2.44721 2.44721  2.4472
## A.het_curt  1 1.39873 1.39873  1.3987
## 
## $FRMA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq  F value
## site        2  80.363  40.182  40.1817
## pesticide   3   0.819   0.273   0.2729
## exclosure   1   3.003   3.003   3.0031
## census      1 177.048 177.048 177.0479
## A.con_curt  1   0.774   0.774   0.7744
## A.het_curt  1   4.268   4.268   4.2678
## 
## $PIKO
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 2.3499  1.1749  1.1749
## pesticide   3 3.0691  1.0230  1.0230
## exclosure   1 2.7894  2.7894  2.7894
## census      1 3.2443  3.2443  3.2443
## A.con_curt  1 1.6280  1.6280  1.6280
## A.het_curt  1 0.1998  0.1998  0.1998
## 
## $QUMO
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 3.7176 1.85881  1.8588
## pesticide   3 3.6884 1.22945  1.2295
## exclosure   1 0.0554 0.05535  0.0554
## census      1 0.0109 0.01088  0.0109
## A.con_curt  1 0.6959 0.69586  0.6959
## A.het_curt  1 1.0059 1.00590  1.0059
## 
## $TIAM
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 12.745   6.372  6.3724
## pesticide   3 29.856   9.952  9.9521
## exclosure   1  6.472   6.472  6.4720
## census      1 83.198  83.198 83.1976
## A.con_curt  1  2.744   2.744  2.7437
## A.het_curt  1  0.525   0.525  0.5248
###----------------- Pooling other species -------------------------###

# sp.rect.sel2 <- c('ACBA',  'ACPS', 'ACTE', 'FRMA', 'PIKO','QUMO', 'TIAM')
pest.t.rect.other_adult <- pest.t.rect_adult[!(pest.t.rect_adult$sp. %in% sp.rect.sel2),]
summary(pest.t.rect.other_adult)
##        X          site      quadrat           sp.     pesticide exclosure
##  Min.   :  21.0   1:44   Min.   :101.0   ACMO   :70   W:39      0:65     
##  1st Qu.: 406.8   2:92   1st Qu.:206.0   ABNE   :45   F:46      1:83     
##  Median : 587.5   3:12   Median :302.5   ACMA   :19   I:37               
##  Mean   : 620.1          Mean   :304.4   MAAM   : 6   C:26               
##  3rd Qu.: 792.5          3rd Qu.:408.0   ULJA   : 3                      
##  Max.   :1424.0          Max.   :510.0   ACTR   : 1                      
##                                          (Other): 4                      
##  growth.form  census      recruits       quad.unique   plot   
##  shrub:  0   16sp:49   Min.   : 1.000   2_0_204:  4   1_0:17  
##  tree :148   17sp:99   1st Qu.: 1.000   2_1_209:  4   1_1:27  
##                        Median : 1.000   1_1_503:  3   2_0:45  
##                        Mean   : 1.743   2_0_206:  3   2_1:47  
##                        3rd Qu.: 2.000   2_0_504:  3   3_0: 3  
##                        Max.   :12.000   2_1_206:  3   3_1: 9  
##                                         (Other):128           
##          quad.sp        A.sum          A.con            A.het     
##  1_0_103_ACMO:  2   Min.   :4106   Min.   :   0.0   Min.   :3560  
##  1_1_503_ACMO:  2   1st Qu.:4949   1st Qu.: 120.3   1st Qu.:4640  
##  2_0_504_ACMO:  2   Median :5397   Median : 268.4   Median :5182  
##  2_1_105_ACMO:  2   Mean   :5620   Mean   : 294.2   Mean   :5326  
##  2_1_209_ULJA:  2   3rd Qu.:6249   3rd Qu.: 394.1   3rd Qu.:5922  
##  2_1_302_ACMO:  2   Max.   :8408   Max.   :1264.4   Max.   :8040  
##  (Other)     :136                                                 
##     con.dens       A.con_curt       A.het_curt      A.sum_curt   
##  Min.   : 0.00   Min.   : 0.000   Min.   :15.27   Min.   :16.01  
##  1st Qu.: 4.75   1st Qu.: 4.937   1st Qu.:16.68   1st Qu.:17.04  
##  Median : 9.00   Median : 6.450   Median :17.31   Median :17.54  
##  Mean   :13.53   Mean   : 5.834   Mean   :17.40   Mean   :17.73  
##  3rd Qu.:20.00   3rd Qu.: 7.331   3rd Qu.:18.09   3rd Qu.:18.42  
##  Max.   :46.00   Max.   :10.813   Max.   :20.03   Max.   :20.33  
##                                                                  
##  con.dens_curt   sp.group
##  Min.   :0.000   0:148   
##  1st Qu.:1.679   1:  0   
##  Median :2.080           
##  Mean   :2.043           
##  3rd Qu.:2.714           
##  Max.   :3.583           
## 
mod.pest.recr.other <- glmer(recruits ~ site + census + pesticide + exclosure + A.con_curt + A.het_curt +
                               (1|quad.unique) + (1|sp.),
                             data=pest.t.rect.other_adult, family=poisson, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0105948 (tol =
## 0.001, component 1)
mod.pest.recr.other.diags <- augment(mod.pest.recr.other)
qplot(data=mod.pest.recr.other.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.pest.recr.other.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.pest.recr.other))

sjp.lmer(mod.pest.recr.other, type='fe')

sjp.lmer(mod.pest.recr.other, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.pest.recr.other, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + census + pesticide + exclosure + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: pest.t.rect.other_adult
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    459.6    495.5   -217.8    435.6      136 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5560 -0.3930 -0.1669  0.1491  4.6519 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  quad.unique (Intercept) 3.572e-05 0.005977
##  sp.         (Intercept) 5.487e-02 0.234235
## Number of obs: 148, groups:  quad.unique, 103; sp., 10
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.832636   0.007448  246.06  < 2e-16 ***
## site1       -0.007413   0.007682   -0.96    0.335    
## site2        0.206950   0.007441   27.81  < 2e-16 ***
## census17sp   0.241130   0.007686   31.37  < 2e-16 ***
## pesticideF  -0.179413   0.007683  -23.35  < 2e-16 ***
## pesticideI   0.063588   0.007440    8.55  < 2e-16 ***
## pesticideC   0.397282   0.007444   53.37  < 2e-16 ***
## exclosure1   0.165368   0.007680   21.53  < 2e-16 ***
## A.con_curt   0.044311   0.007129    6.22 5.11e-10 ***
## A.het_curt  -0.120519   0.005470  -22.03  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.0105948 (tol = 0.001, component 1)

Overall, fungicide higher recruitment, conspecific neigboring trees had a positive effect on recruitment. Tested at species level, fungicide increased recruitment of ACPS, TIAM. Conspecific neighboring trees increased recruitment of ACPS, slightly enhanced TIAM. Heterospecific neighboring trees increased recruitment of FRMA but decreased ACPS. Exclosure enhanced recruitment of ACPS. Pesticide had no effect on recruitment of the rest of species.

Growth analysis

Growth is calculated as the absolute values from the census to the next. SND are used to calculate the growth of seedlings, and outliers are estimated and removed in the final dataset.

## read in data
pestdat.g_adult <- read.csv("data/pestdat.g_adult.csv")
summary(pestdat.g_adult)
##        X              sp.            site          quadrat     
##  Min.   :  170   FRMA   :3144   Min.   :1.000   Min.   :101.0  
##  1st Qu.: 8806   TIAM   :2006   1st Qu.:1.000   1st Qu.:202.0  
##  Median :13889   ACPS   : 571   Median :2.000   Median :304.0  
##  Mean   :13736   EUMA   : 505   Mean   :1.908   Mean   :305.7  
##  3rd Qu.:19216   PHSC   : 484   3rd Qu.:2.000   3rd Qu.:409.0  
##  Max.   :23888   ACBA   : 355   Max.   :3.000   Max.   :510.0  
##                  (Other):2206                                  
##       tag        pesticide   exclosure      growth.form   census    
##  Min.   : 1001   C:1986    Min.   :0.0000   shrub:2452   15fa: 778  
##  1st Qu.: 3018   F:2816    1st Qu.:0.0000   tree :6819   16fa:2214  
##  Median : 6005   I:2177    Median :1.0000                16sp:3088  
##  Mean   : 5701   W:2292    Mean   :0.5113                17sp:3191  
##  3rd Qu.: 8023             3rd Qu.:1.0000                           
##  Max.   :10052             Max.   :1.0000                           
##                                                                     
##     grow.snd          quad.unique    plot              quad.sp    
##  Min.   :-2.211138   2_0_510: 106   1_0:1387   2_1_310_FRMA:  52  
##  1st Qu.:-0.406868   2_1_310:  96   1_1:1703   2_1_103_FRMA:  48  
##  Median : 0.015854   3_0_108:  87   2_0:1967   3_0_108_FRMA:  48  
##  Mean   : 0.002379   2_1_505:  84   2_1:1981   3_1_510_FRMA:  46  
##  3rd Qu.: 0.347986   2_1_103:  79   3_0:1177   2_1_506_FRMA:  44  
##  Max.   : 2.249671   2_1_301:  79   3_1:1056   2_1_304_FRMA:  43  
##                      (Other):8740              (Other)     :8990  
##      A.sum           A.con            A.het          con.dens     
##  Min.   : 3990   Min.   :   0.0   Min.   : 2374   Min.   : 0.000  
##  1st Qu.: 5026   1st Qu.:   0.0   1st Qu.: 4296   1st Qu.: 0.000  
##  Median : 5554   Median : 276.8   Median : 5097   Median : 3.000  
##  Mean   : 5719   Mean   : 593.3   Mean   : 5125   Mean   : 6.032  
##  3rd Qu.: 6213   3rd Qu.: 977.7   3rd Qu.: 5791   3rd Qu.:10.000  
##  Max.   :10411   Max.   :5308.9   Max.   :10411   Max.   :46.000  
##                                                                   
##    A.con_curt       A.het_curt      A.sum_curt    con.dens_curt  
##  Min.   : 0.000   Min.   :13.34   Min.   :15.86   Min.   :0.000  
##  1st Qu.: 0.000   1st Qu.:16.26   1st Qu.:17.13   1st Qu.:0.000  
##  Median : 6.517   Median :17.21   Median :17.71   Median :1.442  
##  Mean   : 5.734   Mean   :17.13   Mean   :17.82   Mean   :1.239  
##  3rd Qu.: 9.925   3rd Qu.:17.96   3rd Qu.:18.38   3rd Qu.:2.154  
##  Max.   :17.445   Max.   :21.84   Max.   :21.84   Max.   :3.583  
##                                                                  
##     season          year      
##  fall  :2992   Min.   :1.000  
##  spring:6279   1st Qu.:1.000  
##                Median :2.000  
##                Mean   :1.583  
##                3rd Qu.:2.000  
##                Max.   :2.000  
## 
###----------------------- Growth models ------------------------------###
pestdat.g_adult$pesticide <- factor(pestdat.g_adult$pesticide, levels=c('W', 'F', 'I', 'C'))
pestdat.g_adult$census <- factor(pestdat.g_adult$census)
pestdat.g_adult$exclosure <- factor(pestdat.g_adult$exclosure)
pestdat.g_adult$site <- factor(pestdat.g_adult$site)
contrasts(pestdat.g_adult$site) <- "contr.sum"

## tree
mod.pest.tree.grow1 <- lmer(grow.snd ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
                         (1|quad.unique) + (1|sp.), 
                       data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow2 <- lmer(grow.snd ~ site + pesticide + exclosure + census + A.con_curt + A.het_curt +
                         (1|quad.unique) + (1|sp.) + (1|tag), 
                       data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow3 <- lmer(grow.snd ~ site + pesticide + exclosure + census +  
                         (1|quad.unique) + (1|sp.), 
                       data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow4 <- lmer(grow.snd ~ site + pesticide + exclosure +  
                              (1|quad.unique) + (1|census) + (1|sp.), 
                            data=subset(pestdat.g_adult, growth.form == 'tree'))

anova(mod.pest.tree.grow1, mod.pest.tree.grow2, mod.pest.tree.grow3, mod.pest.tree.grow4)
## refitting model(s) with ML (instead of REML)
## Data: subset(pestdat.g_adult, growth.form == "tree")
## Models:
## mod.pest.tree.grow4: grow.snd ~ site + pesticide + exclosure + (1 | quad.unique) + 
## mod.pest.tree.grow4:     (1 | census) + (1 | sp.)
## mod.pest.tree.grow3: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) + 
## mod.pest.tree.grow3:     (1 | sp.)
## mod.pest.tree.grow1: grow.snd ~ site + pesticide + exclosure + census + A.con_curt + 
## mod.pest.tree.grow1:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow2: grow.snd ~ site + pesticide + exclosure + census + A.con_curt + 
## mod.pest.tree.grow2:     A.het_curt + (1 | quad.unique) + (1 | sp.) + (1 | tag)
##                     Df   AIC   BIC  logLik deviance  Chisq Chi Df
## mod.pest.tree.grow4 11 13862 13937 -6919.8    13840              
## mod.pest.tree.grow3 13 13843 13932 -6908.5    13817 22.676      2
## mod.pest.tree.grow1 15 13830 13933 -6900.2    13800 16.609      2
## mod.pest.tree.grow2 16 13832 13942 -6900.2    13800  0.000      1
##                     Pr(>Chisq)    
## mod.pest.tree.grow4               
## mod.pest.tree.grow3  1.191e-05 ***
## mod.pest.tree.grow1  0.0002474 ***
## mod.pest.tree.grow2  1.0000000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.pest.tree.grow.diags <- augment(mod.pest.tree.grow1)
qplot(data=mod.pest.tree.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.pest.tree.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.pest.tree.grow1))

summary(mod.pest.tree.grow1, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + pesticide + exclosure + census + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: subset(pestdat.g_adult, growth.form == "tree")
## 
## REML criterion at convergence: 13876.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8929 -0.5715 -0.0315  0.5488  3.4367 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.004908 0.07006 
##  sp.         (Intercept) 0.004189 0.06472 
##  Residual                0.438786 0.66241 
## Number of obs: 6819, groups:  quad.unique, 238; sp., 19
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.165058   0.165519   0.997
## site1        0.022380   0.014858   1.506
## site2        0.020914   0.013392   1.562
## pesticideF   0.004157   0.025929   0.160
## pesticideI  -0.023546   0.027367  -0.860
## pesticideC  -0.045913   0.028072  -1.636
## exclosure1   0.004993   0.019784   0.252
## census16fa   0.130424   0.047378   2.753
## census16sp  -0.347956   0.046609  -7.465
## census17sp  -0.504236   0.046588 -10.823
## A.con_curt  -0.009241   0.003210  -2.879
## A.het_curt   0.011426   0.008806   1.298
###------------------- season and year effect? --------------------------###
unique(pestdat.g_adult$census)
## [1] 15fa 16sp 16fa 17sp
## Levels: 15fa 16fa 16sp 17sp
pestdat.g_adult$season <- ifelse(pestdat.g_adult$census=="16sp" | 
                                   pestdat.g_adult$census=="17sp", "spring","fall")
pestdat.g_adult$year <- ifelse(pestdat.g_adult$census=="15fa" | 
                                   pestdat.g_adult$census=="16sp", "1","2")

str(pestdat.g_adult)
## 'data.frame':    9271 obs. of  23 variables:
##  $ X            : int  170 174 184 248 253 273 280 293 297 298 ...
##  $ sp.          : Factor w/ 38 levels "ABNE","ACBA",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ site         : Factor w/ 3 levels "1","2","3": 2 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "contrasts")= chr "contr.sum"
##  $ quadrat      : int  402 108 309 304 304 209 310 505 503 308 ...
##  $ tag          : int  2003 8005 9007 4002 4011 9008 10007 5002 3004 8002 ...
##  $ pesticide    : Factor w/ 4 levels "W","F","I","C": 1 1 1 3 3 2 2 3 2 4 ...
##  $ exclosure    : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 2 1 1 2 ...
##  $ growth.form  : Factor w/ 2 levels "shrub","tree": 2 2 2 2 2 2 2 2 2 2 ...
##  $ census       : Factor w/ 4 levels "15fa","16fa",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ grow.snd     : num  0.2773 -0.0239 -0.2248 1.9844 0.9802 ...
##  $ quad.unique  : Factor w/ 238 levels "1_0_101","1_0_102",..: 145 8 63 20 20 55 64 36 34 62 ...
##  $ plot         : Factor w/ 6 levels "1_0","1_1","2_0",..: 4 1 2 1 1 2 2 1 1 2 ...
##  $ quad.sp      : Factor w/ 1417 levels "1_0_101_ACKO",..: 997 54 414 135 135 361 424 230 216 405 ...
##  $ A.sum        : num  6306 5299 5626 4975 4975 ...
##  $ A.con        : num  2.89 2.79 1.85 7.86 7.86 ...
##  $ A.het        : num  6303 5296 5624 4967 4967 ...
##  $ con.dens     : int  2 1 1 1 1 1 2 2 1 1 ...
##  $ A.con_curt   : num  1.42 1.41 1.23 1.99 1.99 ...
##  $ A.het_curt   : num  18.5 17.4 17.8 17.1 17.1 ...
##  $ A.sum_curt   : num  18.5 17.4 17.8 17.1 17.1 ...
##  $ con.dens_curt: num  1.26 1 1 1 1 ...
##  $ season       : chr  "fall" "fall" "fall" "fall" ...
##  $ year         : chr  "1" "1" "1" "1" ...
mod.pest.tree.grow5 <- lmer(grow.snd ~ site + pesticide + exclosure + season*year + A.con_curt + A.het_curt +
                              (1|quad.unique) + (1|sp.), 
                            data=subset(pestdat.g_adult, growth.form == 'tree'))
mod.pest.tree.grow6 <- lmer(grow.snd ~ site + pesticide + exclosure + season + year + A.con_curt + A.het_curt +
                              (1|quad.unique) + (1|sp.), 
                            data=subset(pestdat.g_adult, growth.form == 'tree'))

anova(mod.pest.tree.grow1, mod.pest.tree.grow5, mod.pest.tree.grow6)
## refitting model(s) with ML (instead of REML)
## Data: subset(pestdat.g_adult, growth.form == "tree")
## Models:
## mod.pest.tree.grow6: grow.snd ~ site + pesticide + exclosure + season + year + A.con_curt + 
## mod.pest.tree.grow6:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow1: grow.snd ~ site + pesticide + exclosure + census + A.con_curt + 
## mod.pest.tree.grow1:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.pest.tree.grow5: grow.snd ~ site + pesticide + exclosure + season * year + A.con_curt + 
## mod.pest.tree.grow5:     A.het_curt + (1 | quad.unique) + (1 | sp.)
##                     Df   AIC   BIC  logLik deviance  Chisq Chi Df
## mod.pest.tree.grow6 14 13860 13955 -6915.8    13832              
## mod.pest.tree.grow1 15 13830 13933 -6900.2    13800 31.225      1
## mod.pest.tree.grow5 15 13830 13933 -6900.2    13800  0.000      0
##                     Pr(>Chisq)    
## mod.pest.tree.grow6               
## mod.pest.tree.grow1  2.297e-08 ***
## mod.pest.tree.grow5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.pest.tree.grow5, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + season * year + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | sp.)
##    Data: subset(pestdat.g_adult, growth.form == "tree")
## 
## REML criterion at convergence: 13876.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8929 -0.5715 -0.0315  0.5488  3.4367 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.004908 0.07006 
##  sp.         (Intercept) 0.004189 0.06472 
##  Residual                0.438786 0.66241 
## Number of obs: 6819, groups:  quad.unique, 238; sp., 19
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)         0.165058   0.165519   0.997
## site1               0.022380   0.014858   1.506
## site2               0.020914   0.013392   1.562
## pesticideF          0.004157   0.025929   0.160
## pesticideI         -0.023546   0.027367  -0.860
## pesticideC         -0.045913   0.028072  -1.636
## exclosure1          0.004993   0.019784   0.252
## seasonspring       -0.347956   0.046609  -7.465
## year2               0.130424   0.047378   2.753
## A.con_curt         -0.009241   0.003210  -2.879
## A.het_curt          0.011426   0.008806   1.298
## seasonspring:year2 -0.286705   0.051323  -5.586
## shrub
mod.pest.shrub.grow1 <- lmer(grow.snd ~ site + pesticide + exclosure + A.sum_curt + census + 
                         (1|quad.unique) + (1|sp.), 
                       data=subset(pestdat.g_adult, growth.form == 'shrub'))
mod.pest.shrub.grow2 <- lmer(grow.snd ~ site + pesticide + exclosure + census + 
                               (1|quad.unique) + (1|sp.), 
                             data=subset(pestdat.g_adult, growth.form == 'shrub'))
mod.pest.shrub.grow3 <- lmer(grow.snd ~ site + pesticide + exclosure +census + 
                         (1|quad.unique) + (1|sp.) + (1|tag), 
                       data=subset(pestdat.g_adult, growth.form == 'shrub'))

anova(mod.pest.shrub.grow1, mod.pest.shrub.grow2, mod.pest.shrub.grow3)
## refitting model(s) with ML (instead of REML)
## Data: subset(pestdat.g_adult, growth.form == "shrub")
## Models:
## mod.pest.shrub.grow2: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) + 
## mod.pest.shrub.grow2:     (1 | sp.)
## mod.pest.shrub.grow1: grow.snd ~ site + pesticide + exclosure + A.sum_curt + census + 
## mod.pest.shrub.grow1:     (1 | quad.unique) + (1 | sp.)
## mod.pest.shrub.grow3: grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) + 
## mod.pest.shrub.grow3:     (1 | sp.) + (1 | tag)
##                      Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.pest.shrub.grow2 13 5139.8 5215.3 -2556.9   5113.8              
## mod.pest.shrub.grow1 14 5141.7 5222.9 -2556.8   5113.7 0.1453      1
## mod.pest.shrub.grow3 14 5141.8 5223.1 -2556.9   5113.8 0.0000      0
##                      Pr(>Chisq)
## mod.pest.shrub.grow2           
## mod.pest.shrub.grow1      0.703
## mod.pest.shrub.grow3      1.000
mod.pest.shrub.grow.diags <- augment(mod.pest.shrub.grow2)
qplot(data=mod.pest.shrub.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.pest.shrub.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.pest.shrub.grow2))

summary(mod.pest.shrub.grow2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + census + (1 | quad.unique) +  
##     (1 | sp.)
##    Data: subset(pestdat.g_adult, growth.form == "shrub")
## 
## REML criterion at convergence: 5167.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5074 -0.4832 -0.0120  0.5113  3.4326 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  sp.         (Intercept) 0.0000   0.0000  
##  Residual                0.4732   0.6879  
## Number of obs: 2452, groups:  quad.unique, 192; sp., 19
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.168721   0.040591   4.157
## site1       -0.001252   0.019558  -0.064
## site2       -0.007465   0.020397  -0.366
## pesticideF  -0.023955   0.038375  -0.624
## pesticideI   0.019500   0.039790   0.490
## pesticideC  -0.014190   0.039680  -0.358
## exclosure1   0.014623   0.027928   0.524
## census16fa   0.057311   0.040608   1.411
## census16sp  -0.272208   0.040262  -6.761
## census17sp  -0.342546   0.039550  -8.661
##------------------------- individual species -------------------------------##
## And now pulling the species apart. 
## First need to filter down to species that can be analysed - need enough
## reps (lets say 5) and variation in abundance (at least 2x 
sp.grow.counts <- summarise(group_by(pestdat.g_adult, sp.), 
                            n.indiv=length(grow.snd), n.quadrat=length(unique(grow.snd)))

sp.grow.sel <- sp.grow.counts$sp.[sp.grow.counts$n.indiv > 50 & sp.grow.counts$n.quadrat > 10]
sp.grow.counts[sp.grow.counts$sp. %in% sp.grow.sel,]
## # A tibble: 19 × 3
##       sp. n.indiv n.quadrat
##    <fctr>   <int>     <int>
## 1    ACBA     355        29
## 2    ACKO      81        23
## 3    ACMA      78        22
## 4    ACMO     142        17
## 5    ACPS     571        22
## 6    ACTE     190        20
## 7    CEOR     137        31
## 8    DEGL     150        32
## 9    EUAL     231        32
## 10   EUMA     505        36
## 11   FRMA    3144        32
## 12   LOSP      94        26
## 13   PHSC     484        45
## 14   QUMO      80        18
## 15   RIMA      80        26
## 16   RIMX     158        40
## 17   SPCH     206        42
## 18   SYRE     167        28
## 19   TIAM    2006        24
sp.grow.sel2 <- c('ACBA','ACKO','ACMA','ACMO','ACPS','ACTE','CEOR','DEGL','EUAL','EUMA', 
                  'FRMA','LOSP','PHSC','QUMO','RIMX','SPCH','SYRE', 'TIAM')

growmods <- sapply(sp.grow.sel2, function(sp){
  print(sp)
  spdat <- filter(pestdat.g_adult, sp. == sp)
  spdat <- droplevels(spdat)
  mod <- lmer(grow.snd ~ site + pesticide + exclosure + 
                 A.con_curt + A.het_curt + census + (1|quad.unique), 
               data=spdat)
}, simplify=FALSE)
## [1] "ACBA"
## [1] "ACKO"
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## [1] "ACMA"
## [1] "ACMO"
## [1] "ACPS"
## [1] "ACTE"
## [1] "CEOR"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "DEGL"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "EUAL"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "EUMA"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "FRMA"
## [1] "LOSP"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "PHSC"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "QUMO"
## [1] "RIMX"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "SPCH"
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## [1] "SYRE"
## [1] "TIAM"
## diagnostic plots
indmods.grow.resids <- lapply(growmods, augment) 

indmods.grow.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.grow.resids)), dat=indmods.grow.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.grow.resids, aes(x=.fitted, y=.wtres, colour=pesticide)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'

library(lattice)
lapply(growmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACBA
## $ACBA$quad.unique

## 
## 
## $ACKO
## $ACKO$quad.unique

## 
## 
## $ACMA
## $ACMA$quad.unique

## 
## 
## $ACMO
## $ACMO$quad.unique

## 
## 
## $ACPS
## $ACPS$quad.unique

## 
## 
## $ACTE
## $ACTE$quad.unique

## 
## 
## $CEOR
## $CEOR$quad.unique

## 
## 
## $DEGL
## $DEGL$quad.unique

## 
## 
## $EUAL
## $EUAL$quad.unique

## 
## 
## $EUMA
## $EUMA$quad.unique

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## 
## $LOSP
## $LOSP$quad.unique

## 
## 
## $PHSC
## $PHSC$quad.unique

## 
## 
## $QUMO
## $QUMO$quad.unique

## 
## 
## $RIMX
## $RIMX$quad.unique

## 
## 
## $SPCH
## $SPCH$quad.unique

## 
## 
## $SYRE
## $SYRE$quad.unique

## 
## 
## $TIAM
## $TIAM$quad.unique

lapply(growmods, function(x) summary(x, correlation=FALSE))
## $ACBA
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 118.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8440 -0.3841 -0.0551  0.2699  6.0125 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.001069 0.0327  
##  Residual                0.070515 0.2655  
## Number of obs: 355, groups:  quad.unique, 96
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.02986    0.40147  -0.074
## site2        0.01970    0.04262   0.462
## site3       -0.02626    0.05608  -0.468
## pesticideF   0.03613    0.04486   0.805
## pesticideI   0.06464    0.04259   1.518
## pesticideC   0.03307    0.04386   0.754
## exclosure1  -0.03020    0.03619  -0.835
## A.con_curt   0.02986    0.03113   0.959
## A.het_curt   0.01474    0.02201   0.670
## census16fa  -0.05926    0.07281  -0.814
## census16sp  -0.23297    0.07054  -3.303
## census17sp  -0.24286    0.06951  -3.494
## 
## $ACKO
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 123.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2651 -0.5792 -0.0307  0.5207  2.9520 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.008899 0.09433 
##  Residual                0.229533 0.47910 
## Number of obs: 81, groups:  quad.unique, 14
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.17891    1.59398  -0.112
## site2        0.20542    0.32586   0.630
## site3        0.24281    0.46849   0.518
## pesticideF  -0.14537    0.18264  -0.796
## pesticideI  -0.22591    0.19898  -1.135
## pesticideC   0.02469    0.19601   0.126
## A.het_curt   0.01079    0.09073   0.119
## census16fa  -0.15840    0.15464  -1.024
## census16sp  -0.07924    0.14989  -0.529
## census17sp  -0.15407    0.14989  -1.028
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## 
## $ACMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 153.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.13201 -0.52436  0.07435  0.45412  2.60520 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  Residual                0.3677   0.6064  
## Number of obs: 78, groups:  quad.unique, 28
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  2.083091   1.477054   1.410
## site2        0.074927   0.169845   0.441
## site3        0.080634   0.656902   0.123
## pesticideF   0.048753   0.226693   0.215
## pesticideI   0.106156   0.238536   0.445
## pesticideC   0.208450   0.300643   0.693
## exclosure1   0.135532   0.198238   0.684
## A.con_curt   0.006269   0.071242   0.088
## A.het_curt  -0.144097   0.087893  -1.640
## census16fa   0.600047   0.230200   2.607
## census16sp   0.102611   0.221310   0.464
## census17sp   0.005543   0.226436   0.024
## 
## $ACMO
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 289.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2511 -0.5725  0.0954  0.6358  2.9564 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  Residual                0.3979   0.6308  
## Number of obs: 142, groups:  quad.unique, 64
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.61605    1.28091   0.481
## site2        0.20873    0.14123   1.478
## site3       -0.13429    0.18291  -0.734
## pesticideF  -0.16813    0.14258  -1.179
## pesticideI  -0.17208    0.16448  -1.046
## pesticideC  -0.14416    0.19728  -0.731
## exclosure1  -0.09585    0.11379  -0.842
## A.con_curt   0.05493    0.06117   0.898
## A.het_curt  -0.03423    0.06156  -0.556
## census16fa  -0.20644    0.39926  -0.517
## census16sp  -0.34968    0.39367  -0.888
## census17sp  -0.41797    0.38980  -1.072
## 
## $ACPS
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 935.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2362 -0.5511 -0.0253  0.6237  4.0505 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 1.180e-15 3.435e-08
##  Residual                2.818e-01 5.308e-01
## Number of obs: 571, groups:  quad.unique, 146
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.65122    0.47133   1.382
## site2        0.08588    0.05462   1.572
## site3       -0.20869    0.06869  -3.038
## pesticideF   0.03595    0.05972   0.602
## pesticideI  -0.02598    0.06450  -0.403
## pesticideC   0.01486    0.06975   0.213
## exclosure1  -0.00874    0.05560  -0.157
## A.con_curt   0.02879    0.02960   0.972
## A.het_curt  -0.02989    0.02548  -1.173
## census16fa   0.25725    0.15374   1.673
## census16sp  -0.34192    0.15282  -2.237
## census17sp  -0.36838    0.14550  -2.532
## 
## $ACTE
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 205.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0577 -0.5037 -0.0451  0.4431  4.2983 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  Residual                0.1467   0.3831  
## Number of obs: 190, groups:  quad.unique, 48
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) -1.009530   0.688932  -1.465
## site2        0.002721   0.069928   0.039
## pesticideF  -0.032059   0.084703  -0.378
## pesticideI   0.068001   0.075892   0.896
## pesticideC  -0.021987   0.084745  -0.260
## exclosure1   0.212109   0.088453   2.398
## A.con_curt   0.070913   0.040015   1.772
## A.het_curt   0.033758   0.037922   0.890
## census16fa   0.266343   0.154524   1.724
## census16sp   0.055113   0.151960   0.363
## census17sp  -0.024637   0.150177  -0.164
## 
## $CEOR
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 304.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9594 -0.3373  0.0166  0.4735  2.6755 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 2.695e-15 5.191e-08
##  Residual                4.981e-01 7.058e-01
## Number of obs: 137, groups:  quad.unique, 22
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  1.382241   1.500046   0.922
## site2        0.072388   0.156682   0.462
## pesticideF   0.132996   0.209171   0.636
## pesticideI   0.278383   0.174202   1.598
## pesticideC  -0.012866   0.192453  -0.067
## exclosure1   0.243905   0.440942   0.553
## A.het_curt  -0.074669   0.087294  -0.855
## census16fa   0.006404   0.181177   0.035
## census16sp  -0.265446   0.178624  -1.486
## census17sp  -0.473348   0.175663  -2.695
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $DEGL
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 339.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3164 -0.4745  0.0839  0.5065  2.9166 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.005447 0.0738  
##  Residual                0.509337 0.7137  
## Number of obs: 150, groups:  quad.unique, 34
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.81802    1.09937  -0.744
## site2        0.10103    0.18156   0.556
## site3        0.02975    0.17806   0.167
## pesticideF   0.03134    0.18186   0.172
## pesticideI  -0.12858    0.17152  -0.750
## pesticideC  -0.25382    0.19298  -1.315
## exclosure1  -0.13297    0.17756  -0.749
## A.het_curt   0.07245    0.06605   1.097
## census16fa   0.01018    0.16870   0.060
## census16sp  -0.61525    0.16735  -3.676
## census17sp  -0.59917    0.16564  -3.617
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $EUAL
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 540.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2848 -0.5162 -0.0112  0.6156  2.4034 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.00    
##  Residual                0.5624   0.75    
## Number of obs: 231, groups:  quad.unique, 46
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  2.078004   1.226454   1.694
## site2       -0.039674   0.118845  -0.334
## site3       -0.179866   0.159746  -1.126
## pesticideF   0.001971   0.149129   0.013
## pesticideI   0.157589   0.160361   0.983
## pesticideC   0.007698   0.147744   0.052
## exclosure1   0.034351   0.116000   0.296
## A.het_curt  -0.102106   0.069144  -1.477
## census16fa   0.255233   0.144946   1.761
## census16sp  -0.513113   0.145435  -3.528
## census17sp  -0.572619   0.140878  -4.065
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $EUMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 1065.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.11938 -0.48812 -0.01841  0.60278  2.85305 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.000    0.0000  
##  Residual                0.456    0.6753  
## Number of obs: 505, groups:  quad.unique, 81
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.38955    0.51048   0.763
## site2        0.01009    0.07253   0.139
## site3        0.22049    0.11242   1.961
## pesticideF   0.06040    0.08405   0.719
## pesticideI   0.11650    0.08924   1.306
## pesticideC   0.15739    0.08716   1.806
## exclosure1  -0.05267    0.06629  -0.795
## A.het_curt  -0.02556    0.02988  -0.856
## census16fa   0.26496    0.08803   3.010
## census16sp  -0.08835    0.08781  -1.006
## census17sp  -0.08136    0.08592  -0.947
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $FRMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 6175.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1668 -0.6048  0.0217  0.6321  3.8548 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 1.027e-15 3.204e-08
##  Residual                4.097e-01 6.401e-01
## Number of obs: 3144, groups:  quad.unique, 224
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.376243   0.184366   2.041
## site2       -0.080184   0.033388  -2.402
## site3       -0.086539   0.035897  -2.411
## pesticideF   0.015605   0.031668   0.493
## pesticideI  -0.034717   0.032590  -1.065
## pesticideC  -0.060307   0.032638  -1.848
## exclosure1  -0.026709   0.023596  -1.132
## A.con_curt  -0.009061   0.003407  -2.659
## A.het_curt   0.003696   0.009837   0.376
## census16fa   0.210417   0.056090   3.751
## census16sp  -0.372730   0.055211  -6.751
## census17sp  -0.669800   0.056171 -11.924
## 
## $LOSP
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 205.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2979 -0.5286  0.0695  0.6064  2.1696 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  Residual                0.4736   0.6882  
## Number of obs: 94, groups:  quad.unique, 22
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.90763    1.93517  -0.469
## site2       -0.23221    0.22873  -1.015
## site3       -0.10166    0.19955  -0.509
## pesticideF   0.20728    0.21811   0.950
## pesticideI  -0.10693    0.22169  -0.482
## pesticideC  -0.04189    0.20987  -0.200
## exclosure1  -0.33438    0.17138  -1.951
## A.het_curt   0.07335    0.10588   0.693
## census16fa  -0.07167    0.20664  -0.347
## census16sp  -0.70535    0.20854  -3.382
## census17sp  -0.39325    0.20287  -1.938
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $PHSC
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 1044.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3339 -0.5485 -0.0017  0.5342  3.0945 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  Residual                0.4778   0.6912  
## Number of obs: 484, groups:  quad.unique, 84
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  0.5919165  0.7138356   0.829
## site2       -0.0057196  0.0709174  -0.081
## site3       -0.1024200  0.0975137  -1.050
## pesticideF  -0.0002306  0.0927566  -0.002
## pesticideI   0.0176862  0.1070802   0.165
## pesticideC   0.0467419  0.0979844   0.477
## exclosure1   0.1416345  0.0641147   2.209
## A.het_curt  -0.0209633  0.0413780  -0.507
## census16fa  -0.2154622  0.0910099  -2.367
## census16sp  -0.3959989  0.0893055  -4.434
## census17sp  -0.4719985  0.0885211  -5.332
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $QUMO
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 147.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75492 -0.57954 -0.08754  0.45398  2.25866 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1431   0.3782  
##  Residual                0.2485   0.4985  
## Number of obs: 80, groups:  quad.unique, 38
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.57531    2.62537   0.981
## site2       -0.20281    0.45852  -0.442
## site3        0.05050    0.28576   0.177
## pesticideF  -0.04227    0.24932  -0.170
## pesticideI   0.12223    0.25590   0.478
## pesticideC  -0.12616    0.32078  -0.393
## exclosure1   0.06607    0.32081   0.206
## A.con_curt  -0.08607    0.08279  -1.040
## A.het_curt  -0.08855    0.13988  -0.633
## census16fa  -0.34245    0.58694  -0.584
## census16sp  -0.66438    0.58727  -1.131
## census17sp  -0.39055    0.57485  -0.679
## 
## $RIMX
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 294.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7970 -0.3281  0.0329  0.4977  2.6292 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.000    0.0000  
##  Residual                0.334    0.5779  
## Number of obs: 158, groups:  quad.unique, 35
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -1.35783    1.11062  -1.223
## site2        0.04218    0.10337   0.408
## site3        0.07470    0.14476   0.516
## pesticideF   0.03057    0.11456   0.267
## pesticideI   0.01713    0.12897   0.133
## pesticideC  -0.04483    0.18020  -0.249
## exclosure1  -0.08953    0.09861  -0.908
## A.het_curt   0.08567    0.06278   1.365
## census16fa  -0.02938    0.13604  -0.216
## census16sp  -0.04387    0.13508  -0.325
## census17sp  -0.25488    0.13201  -1.931
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $SPCH
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 405.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1798 -0.5181 -0.0613  0.5629  2.8514 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 1.092e-18 1.045e-09
##  Residual                3.786e-01 6.153e-01
## Number of obs: 206, groups:  quad.unique, 38
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.315818   1.165221   0.271
## site2       -0.078148   0.133826  -0.584
## site3       -0.004012   0.112478  -0.036
## pesticideF  -0.119994   0.141277  -0.849
## pesticideI  -0.258568   0.117254  -2.205
## pesticideC  -0.085017   0.127398  -0.667
## exclosure1  -0.180097   0.100746  -1.788
## A.het_curt   0.009809   0.064988   0.151
## census16fa  -0.205046   0.125477  -1.634
## census16sp  -0.114337   0.123404  -0.927
## census17sp  -0.460315   0.124963  -3.684
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 
## $SYRE
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 370.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4802 -0.4671  0.0078  0.4670  2.8430 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 3.554e-16 1.885e-08
##  Residual                4.853e-01 6.966e-01
## Number of obs: 167, groups:  quad.unique, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -1.00354    1.21507  -0.826
## site2        0.20770    0.17993   1.154
## site3        0.17573    0.20076   0.875
## pesticideF  -0.01741    0.18481  -0.094
## pesticideI   0.30822    0.17605   1.751
## pesticideC  -0.17026    0.18533  -0.919
## exclosure1  -0.08093    0.16180  -0.500
## A.con_curt   0.04516    0.08304   0.544
## A.het_curt   0.05423    0.06468   0.838
## census16fa   0.17167    0.15853   1.083
## census16sp  -0.28124    0.15766  -1.784
## census17sp  -0.28498    0.15404  -1.850
## 
## $TIAM
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique)
##    Data: spdat
## 
## REML criterion at convergence: 4630.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01406 -0.58696 -0.06046  0.55875  2.82722 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.01234  0.1111  
##  Residual                0.56533  0.7519  
## Number of obs: 2006, groups:  quad.unique, 225
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.642327   0.434866   1.477
## site2        0.011795   0.051174   0.230
## site3       -0.080662   0.049428  -1.632
## pesticideF   0.029144   0.050777   0.574
## pesticideI  -0.009912   0.056936  -0.174
## pesticideC   0.005283   0.059945   0.088
## exclosure1   0.015908   0.040225   0.396
## A.con_curt  -0.040216   0.014855  -2.707
## A.het_curt   0.024414   0.018220   1.340
## census16fa  -0.292464   0.244175  -1.198
## census16sp  -0.630650   0.242611  -2.599
## census17sp  -0.726529   0.243035  -2.989
lapply(growmods, function(x) anova(x))
## $ACBA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 0.19044 0.09522  1.3503
## pesticide   3 0.10601 0.03534  0.5011
## exclosure   1 0.04581 0.04581  0.6497
## A.con_curt  1 0.04287 0.04287  0.6080
## A.het_curt  1 0.02771 0.02771  0.3929
## census      3 2.44742 0.81581 11.5692
## 
## $ACKO
## Analysis of Variance Table
##            Df  Sum Sq  Mean Sq F value
## site        2 0.13388 0.066938  0.2916
## pesticide   3 0.56789 0.189298  0.8247
## A.het_curt  1 0.00612 0.006123  0.0267
## census      3 0.33188 0.110625  0.4820
## 
## $ACMA
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 0.1982 0.09908  0.2695
## pesticide   3 0.0500 0.01667  0.0453
## exclosure   1 0.1599 0.15986  0.4348
## A.con_curt  1 0.0155 0.01546  0.0421
## A.het_curt  1 0.6869 0.68691  1.8683
## census      3 4.4915 1.49716  4.0720
## 
## $ACMO
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 1.62732 0.81366  2.0450
## pesticide   3 0.75575 0.25192  0.6332
## exclosure   1 0.16894 0.16894  0.4246
## A.con_curt  1 0.26634 0.26634  0.6694
## A.het_curt  1 0.06992 0.06992  0.1757
## census      3 1.22562 0.40854  1.0268
## 
## $ACPS
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  6.1451  3.0726 10.9050
## pesticide   3  0.6778  0.2259  0.8019
## exclosure   1  0.2705  0.2705  0.9600
## A.con_curt  1  0.2524  0.2524  0.8957
## A.het_curt  1  0.3625  0.3625  1.2865
## census      3 27.7756  9.2585 32.8599
## 
## $ACTE
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        1 0.42257 0.42257  2.8798
## pesticide   3 0.16310 0.05437  0.3705
## exclosure   1 0.53553 0.53553  3.6496
## A.con_curt  1 0.65353 0.65353  4.4538
## A.het_curt  1 0.10551 0.10551  0.7190
## census      3 2.48124 0.82708  5.6365
## 
## $CEOR
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        1 0.0036 0.00358  0.0072
## pesticide   3 1.7231 0.57435  1.1531
## exclosure   1 0.0012 0.00118  0.0024
## A.het_curt  1 0.3364 0.33640  0.6754
## census      3 5.6007 1.86689  3.7480
## 
## $DEGL
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  0.8091  0.4046  0.7943
## pesticide   3  0.9130  0.3043  0.5975
## exclosure   1  0.1115  0.1115  0.2189
## A.het_curt  1  0.6542  0.6542  1.2843
## census      3 13.9627  4.6542  9.1378
## 
## $EUAL
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  1.0750  0.5375  0.9556
## pesticide   3  0.7356  0.2452  0.4360
## exclosure   1  0.0588  0.0588  0.1045
## A.het_curt  1  1.3798  1.3798  2.4533
## census      3 28.5855  9.5285 16.9417
## 
## $EUMA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  1.3905  0.6953  1.5248
## pesticide   3  1.3571  0.4524  0.9921
## exclosure   1  0.4103  0.4103  0.8999
## A.het_curt  1  0.2873  0.2873  0.6302
## census      3 10.5258  3.5086  7.6946
## 
## $FRMA
## Analysis of Variance Table
##            Df Sum Sq Mean Sq  F value
## site        2   4.10   2.050   5.0026
## pesticide   3   2.73   0.909   2.2180
## exclosure   1   0.79   0.787   1.9208
## A.con_curt  1   4.50   4.497  10.9773
## A.het_curt  1   0.08   0.078   0.1901
## census      3 378.44 126.147 307.9073
## 
## $LOSP
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 0.2139 0.10696  0.2259
## pesticide   3 0.4119 0.13731  0.2900
## exclosure   1 1.7564 1.75638  3.7089
## A.het_curt  1 0.2958 0.29584  0.6247
## census      3 7.1127 2.37092  5.0066
## 
## $PHSC
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2  0.6653  0.3327  0.6963
## pesticide   3  0.1744  0.0581  0.1217
## exclosure   1  2.2935  2.2935  4.8003
## A.het_curt  1  0.0692  0.0692  0.1447
## census      3 15.9818  5.3273 11.1500
## 
## $QUMO
## Analysis of Variance Table
##            Df  Sum Sq  Mean Sq F value
## site        2 0.01168 0.005840  0.0235
## pesticide   3 0.17480 0.058267  0.2345
## exclosure   1 0.00529 0.005290  0.0213
## A.con_curt  1 0.14338 0.143378  0.5769
## A.het_curt  1 0.09580 0.095799  0.3855
## census      3 0.54975 0.183250  0.7374
## 
## $RIMX
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 0.33322 0.16661  0.4989
## pesticide   3 0.00321 0.00107  0.0032
## exclosure   1 0.18359 0.18359  0.5497
## A.het_curt  1 0.72180 0.72180  2.1612
## census      3 1.71309 0.57103  1.7098
## 
## $SPCH
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 0.0178 0.00891  0.0235
## pesticide   3 2.0533 0.68443  1.8078
## exclosure   1 1.1318 1.13185  2.9896
## A.het_curt  1 0.0083 0.00828  0.0219
## census      3 5.7822 1.92740  5.0910
## 
## $SYRE
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 0.0373 0.01866  0.0384
## pesticide   3 5.6666 1.88888  3.8922
## exclosure   1 0.3082 0.30825  0.6352
## A.con_curt  1 0.0592 0.05924  0.1221
## A.het_curt  1 0.3271 0.32714  0.6741
## census      3 6.3363 2.11210  4.3522
## 
## $TIAM
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2  1.060  0.5301  0.9377
## pesticide   3  0.502  0.1674  0.2961
## exclosure   1  0.218  0.2176  0.3849
## A.con_curt  1  7.343  7.3426 12.9882
## A.het_curt  1  1.027  1.0271  1.8168
## census      3 51.364 17.1213 30.2855
##------------------------ Other species ------------------------##
pestdat.g.other_adult <- pestdat.g_adult[!(pestdat.g_adult$sp. %in% sp.grow.sel2),]
mod.pest.other.grow <- lmer(grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt + census + 
                              (1|quad.unique) + (1|sp.), data=pestdat.g.other_adult)

mod.pest.other.grow.diags <- augment(mod.pest.other.grow)
qplot(data=mod.pest.other.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.pest.other.grow.diags, x= pesticide, y=.wtres, geom="boxplot", colour=exclosure) 

qqPlot(resid(mod.pest.other.grow))

summary(mod.pest.other.grow, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## grow.snd ~ site + pesticide + exclosure + A.con_curt + A.het_curt +  
##     census + (1 | quad.unique) + (1 | sp.)
##    Data: pestdat.g.other_adult
## 
## REML criterion at convergence: 1202.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8924 -0.5571 -0.0478  0.5790  3.2409 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  quad.unique (Intercept) 0.0007595 0.02756 
##  sp.         (Intercept) 0.0000000 0.00000 
##  Residual                0.6282488 0.79262 
## Number of obs: 492, groups:  quad.unique, 109; sp., 20
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.935784   0.585889   1.597
## site1       -0.004980   0.057059  -0.087
## site2       -0.054944   0.051073  -1.076
## pesticideF  -0.098299   0.107024  -0.918
## pesticideI   0.083325   0.111394   0.748
## pesticideC   0.043735   0.109637   0.399
## exclosure1   0.348849   0.079076   4.412
## A.con_curt   0.004207   0.010270   0.410
## A.het_curt  -0.053095   0.032825  -1.618
## census16fa   0.141441   0.129259   1.094
## census16sp  -0.302241   0.127118  -2.378
## census17sp  -0.331094   0.121582  -2.723

Pesticide had no effect on growth of seedlings. Results are consistent both at community- and species-level. Exclosure increased growth of the species other than the common species.

Diversity analysis

First, trees and shrubs are pooled and estimated the Shannon’ H and Inversimpson indexes. Next, trees and shrubs are estimated seperately.

library(vegan)
## Loading required package: permute
## This is vegan 2.4-4
##----------------------- Overall --------------------------##

divers_woodpest <- summarise(
  group_by(pestdat.ag_adult, site, plot, quadrat, quad.unique, pesticide, exclosure, census), 
  shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))

adult.ngbr.inf <- read.csv("data/adult.ngbr.inf.csv")

divers_woodpest$A.sum <- adult.ngbr.inf$A.sum[match(divers_woodpest$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_woodpest)
##  site     plot        quadrat       quad.unique  pesticide exclosure
##  1:317   1_0:159   Min.   :101.0   1_0_101:  4   W:234     0:474    
##  2:318   1_1:158   1st Qu.:202.0   1_0_102:  4   F:237     1:463    
##  3:302   2_0:159   Median :304.0   1_0_103:  4   I:230              
##          2_1:159   Mean   :304.7   1_0_104:  4   C:236              
##          3_0:156   3rd Qu.:409.0   1_0_105:  4                      
##          3_1:146   Max.   :510.0   1_0_106:  4                      
##                                    (Other):913                      
##   census       shannon          simpson         species      
##  16fa:238   Min.   :0.0000   Min.   :1.000   Min.   : 0.000  
##  16sp:223   1st Qu.:0.7356   1st Qu.:2.000   1st Qu.: 3.000  
##  17fa:238   Median :1.2033   Median :2.941   Median : 4.000  
##  17sp:238   Mean   :1.1393   Mean   :  Inf   Mean   : 4.291  
##             3rd Qu.:1.5591   3rd Qu.:4.167   3rd Qu.: 6.000  
##             Max.   :2.2430   Max.   :  Inf   Max.   :11.000  
##                                                              
##      A.sum      
##  Min.   : 3990  
##  1st Qu.: 4975  
##  Median : 5463  
##  Mean   : 5618  
##  3rd Qu.: 6071  
##  Max.   :10411  
## 
divers_woodpest$A.sum_curt <- divers_woodpest$A.sum^(1/3)


ggplot(divers_woodpest, aes(x=pesticide, y=shannon, colour = exclosure)) +  geom_boxplot()

divers_woodpest$census <- factor(divers_woodpest$census, levels=c('16sp', '16fa', '17sp', '17fa'))

summarise(group_by(divers_woodpest, pesticide, exclosure, census), meanH=mean(shannon), 
          seH=sd(shannon)/sqrt(length(shannon))) %>%
  ggplot(aes(x=pesticide, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census, 
             shape=as.factor(exclosure))) + geom_pointrange() + coord_trans(y="exp") +
  labs(x="pesticide", y="Effective number of species")

### shannon diversity
mod.wood.shannon1 <- lmer(shannon ~ site + pesticide*exclosure + census+ (1|quad.unique), data=divers_woodpest)
mod.wood.shannon2 <- lmer(shannon ~ site + pesticide + exclosure + census + (1|quad.unique), data=divers_woodpest)
mod.wood.shannon3 <- lmer(shannon ~ site + pesticide + exclosure + census + A.sum_curt +
                            (1|quad.unique), data=divers_woodpest)
mod.wood.shannon4 <- lmer(shannon ~ pesticide + exclosure + census + (1|plot/quadrat), data=divers_woodpest)

anova(mod.wood.shannon1, mod.wood.shannon2, mod.wood.shannon3, mod.wood.shannon4)
## refitting model(s) with ML (instead of REML)
## Data: divers_woodpest
## Models:
## mod.wood.shannon4: shannon ~ pesticide + exclosure + census + (1 | plot/quadrat)
## mod.wood.shannon2: shannon ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.shannon3: shannon ~ site + pesticide + exclosure + census + A.sum_curt + 
## mod.wood.shannon3:     (1 | quad.unique)
## mod.wood.shannon1: shannon ~ site + pesticide * exclosure + census + (1 | quad.unique)
##                   Df    AIC    BIC  logLik deviance   Chisq Chi Df
## mod.wood.shannon4 11 518.41 571.67 -248.20   496.41               
## mod.wood.shannon2 12 500.14 558.25 -238.07   476.14 20.2701      1
## mod.wood.shannon3 13 500.94 563.89 -237.47   474.94  1.1952      1
## mod.wood.shannon1 15 503.22 575.86 -236.61   473.22  1.7177      2
##                   Pr(>Chisq)    
## mod.wood.shannon4               
## mod.wood.shannon2  6.724e-06 ***
## mod.wood.shannon3     0.2743    
## mod.wood.shannon1     0.4237    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod.wood.shannon2, type=c('p', 'smooth'))

library(sjPlot)
sjp.lmer(mod.wood.shannon2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wood.shannon2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.wood.shannon2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## shannon ~ site + pesticide + exclosure + census + (1 | quad.unique)
##    Data: divers_woodpest
## 
## REML criterion at convergence: 525.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7845 -0.5281  0.0199  0.5461  3.1071 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.12917  0.3594  
##  Residual                0.05428  0.2330  
## Number of obs: 937, groups:  quad.unique, 239
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.74860    0.05639   13.28
## site1        0.24868    0.03456    7.20
## site2        0.06529    0.03456    1.89
## pesticideF   0.06764    0.06906    0.98
## pesticideI  -0.00947    0.06942   -0.14
## pesticideC  -0.05423    0.06907   -0.79
## exclosure1  -0.13516    0.04895   -2.76
## census16fa   0.53316    0.02184   24.41
## census17sp   0.49461    0.02184   22.64
## census17fa   0.72701    0.02184   33.28
##----------------------- More diversity indices -------------------------##

### Invsimpson
range(divers_woodpest$simpson)
## [1]   1 Inf
divers_woodpest1 <- subset(divers_woodpest, simpson != Inf)

mod.wood.simpson1 <- lmer(simpson ~ site + pesticide*exclosure + census + (1|quad.unique), data=divers_woodpest1)
mod.wood.simpson2 <- lmer(simpson ~ site + pesticide + exclosure + census + (1|quad.unique), data=divers_woodpest1)
mod.wood.simpson3 <- lmer(simpson ~ site + pesticide + exclosure + census + A.sum_curt +
                             (1|quad.unique), data=divers_woodpest1)

anova(mod.wood.simpson1, mod.wood.simpson2, mod.wood.simpson3)
## refitting model(s) with ML (instead of REML)
## Data: divers_woodpest1
## Models:
## mod.wood.simpson2: simpson ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.simpson3: simpson ~ site + pesticide + exclosure + census + A.sum_curt + 
## mod.wood.simpson3:     (1 | quad.unique)
## mod.wood.simpson1: simpson ~ site + pesticide * exclosure + census + (1 | quad.unique)
##                   Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.wood.simpson2 12 2610.7 2668.7 -1293.4   2586.7              
## mod.wood.simpson3 13 2612.4 2675.3 -1293.2   2586.4 0.2661      1
## mod.wood.simpson1 15 2615.6 2688.1 -1292.8   2585.6 0.8460      2
##                   Pr(>Chisq)
## mod.wood.simpson2           
## mod.wood.simpson3     0.6059
## mod.wood.simpson1     0.6551
plot(mod.wood.simpson2, type=c('p', 'smooth'))

sjp.lmer(mod.wood.simpson2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wood.simpson2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.wood.simpson2))

summary(mod.wood.simpson2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## simpson ~ site + pesticide + exclosure + census + (1 | quad.unique)
##    Data: divers_woodpest1
## 
## REML criterion at convergence: 2614.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9611 -0.5662 -0.0780  0.5089  4.1866 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.9858   0.9929  
##  Residual                0.5688   0.7542  
## Number of obs: 928, groups:  quad.unique, 238
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  2.324130   0.160422  14.488
## site1        0.759747   0.097315   7.807
## site2        0.004584   0.097315   0.047
## pesticideF   0.201645   0.194368   1.037
## pesticideI   0.073678   0.196057   0.376
## pesticideC  -0.090029   0.194237  -0.463
## exclosure1  -0.394519   0.138010  -2.859
## census16fa   0.957426   0.071367  13.415
## census17sp   1.012545   0.071462  14.169
## census17fa   1.750211   0.071409  24.510
### Pielou's evenness (J)
table(divers_woodpest$species)
## 
##   0   1   2   3   4   5   6   7   8   9  10  11 
##   9  85 118 181 139 127 116  72  56  24   8   2
# need to remove those with less than 1 species, do not test this index in trees and shrubs seperately.
divers_woodpest2 <- subset(divers_woodpest, species > 1)
divers_woodpest2$evenness <- divers_woodpest2$shannon/log(divers_woodpest2$species)

mod.wood.evenness1 <- lmer(evenness ~ site +  pesticide*exclosure + census +(1|quad.unique), data=divers_woodpest2)
mod.wood.evenness2 <- lmer(evenness ~ site + pesticide + exclosure + census +(1|quad.unique), data=divers_woodpest2)
mod.wood.evenness3 <- lmer(evenness ~ site + pesticide + exclosure + census +A.sum_curt +
                            (1|quad.unique), data=divers_woodpest2)
anova(mod.wood.evenness1, mod.wood.evenness2, mod.wood.evenness3)
## refitting model(s) with ML (instead of REML)
## Data: divers_woodpest2
## Models:
## mod.wood.evenness2: evenness ~ site + pesticide + exclosure + census + (1 | quad.unique)
## mod.wood.evenness3: evenness ~ site + pesticide + exclosure + census + A.sum_curt + 
## mod.wood.evenness3:     (1 | quad.unique)
## mod.wood.evenness1: evenness ~ site + pesticide * exclosure + census + (1 | quad.unique)
##                    Df     AIC     BIC logLik deviance  Chisq Chi Df
## mod.wood.evenness2 12 -1730.1 -1673.2 877.03  -1754.1              
## mod.wood.evenness3 13 -1734.4 -1672.8 880.22  -1760.4 6.3676      1
## mod.wood.evenness1 15 -1726.3 -1655.3 878.17  -1756.3 0.0000      2
##                    Pr(>Chisq)  
## mod.wood.evenness2             
## mod.wood.evenness3    0.01162 *
## mod.wood.evenness1    1.00000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod.wood.evenness2, type=c('p', 'smooth'))

library(sjPlot)
sjp.lmer(mod.wood.evenness2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wood.evenness2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.wood.evenness2))

summary(mod.wood.evenness2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## evenness ~ site + pesticide + exclosure + census + (1 | quad.unique)
##    Data: divers_woodpest2
## 
## REML criterion at convergence: -1676.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4809 -0.5011  0.1263  0.5480  3.3404 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.005098 0.07140 
##  Residual                0.004794 0.06924 
## Number of obs: 843, groups:  quad.unique, 234
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.957464   0.012776   74.94
## site1        0.035170   0.007385    4.76
## site2       -0.037207   0.007385   -5.04
## pesticideF  -0.013875   0.014724   -0.94
## pesticideI  -0.001292   0.014897   -0.09
## pesticideC   0.001681   0.014863    0.11
## exclosure1  -0.037541   0.010535   -3.56
## census16fa  -0.099753   0.007385  -13.51
## census17sp  -0.067521   0.007400   -9.12
## census17fa  -0.058929   0.007393   -7.97
summary(mod.wood.evenness3, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: evenness ~ site + pesticide + exclosure + census + A.sum_curt +  
##     (1 | quad.unique)
##    Data: divers_woodpest2
## 
## REML criterion at convergence: -1674
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5331 -0.5062  0.1335  0.5453  3.3869 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.004952 0.07037 
##  Residual                0.004794 0.06924 
## Number of obs: 843, groups:  quad.unique, 234
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  1.200959   0.098399  12.205
## site1        0.031432   0.007451   4.218
## site2       -0.034372   0.007389  -4.652
## pesticideF  -0.013587   0.014556  -0.933
## pesticideI  -0.000756   0.014728  -0.051
## pesticideC   0.004779   0.014746   0.324
## exclosure1  -0.030991   0.010736  -2.887
## census16fa  -0.099919   0.007384 -13.531
## census17sp  -0.067687   0.007400  -9.147
## census17fa  -0.059049   0.007392  -7.988
## A.sum_curt  -0.013966   0.005598  -2.495
###-------------------- Trees, shrubs -------------------------###
## Split into two growth forms: trees and shrubs
divers_tspest <- summarise(
  group_by(pestdat.ag_adult, site, plot, quadrat, quad.unique, growth.form, pesticide, exclosure, census), 
  shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))

divers_tspest$A.sum <- adult.ngbr.inf$A.sum[match(divers_tspest$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_tspest)
##  site     plot        quadrat       quad.unique   growth.form pesticide
##  1:580   1_0:293   Min.   :101.0   1_0_101:   8   shrub:753   W:400    
##  2:572   1_1:287   1st Qu.:202.0   1_0_103:   8   tree :853   F:407    
##  3:454   2_0:300   Median :304.0   1_0_105:   8               I:394    
##          2_1:272   Mean   :302.9   1_0_107:   8               C:405    
##          3_0:229   3rd Qu.:408.0   1_0_108:   8                        
##          3_1:225   Max.   :510.0   1_0_201:   8                        
##                                    (Other):1558                        
##  exclosure  census       shannon          simpson         species     
##  0:822     16fa:426   Min.   :0.0000   Min.   :1.000   Min.   :0.000  
##  1:784     16sp:327   1st Qu.:0.0000   1st Qu.:1.324   1st Qu.:1.000  
##            17fa:428   Median :0.6931   Median :2.000   Median :2.000  
##            17sp:425   Mean   :0.6933   Mean   :  Inf   Mean   :2.504  
##                       3rd Qu.:1.0735   3rd Qu.:2.830   3rd Qu.:3.000  
##                       Max.   :1.9560   Max.   :  Inf   Max.   :9.000  
##                                                                       
##      A.sum      
##  Min.   : 3990  
##  1st Qu.: 4991  
##  Median : 5455  
##  Mean   : 5630  
##  3rd Qu.: 6045  
##  Max.   :10411  
## 
divers_tspest$A.sum_curt <- divers_tspest$A.sum^(1/3)


### Trees
# shannon diversity
mod.tree.shannon1 <- lmer(shannon ~ site + census + pesticide*exclosure + (1|quad.unique), 
                          data=subset(divers_tspest, growth.form=='tree'))
mod.tree.shannon2 <- lmer(shannon ~ site + census + pesticide + exclosure + (1|quad.unique),
                          data=subset(divers_tspest, growth.form=='tree'))
mod.tree.shannon3 <- lmer(shannon ~ site + census + pesticide + exclosure + A.sum_curt +
                            (1|quad.unique), data=subset(divers_tspest, growth.form=='tree'))
anova(mod.tree.shannon1, mod.tree.shannon2, mod.tree.shannon3)
## refitting model(s) with ML (instead of REML)
## Data: subset(divers_tspest, growth.form == "tree")
## Models:
## mod.tree.shannon2: shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## mod.tree.shannon3: shannon ~ site + census + pesticide + exclosure + A.sum_curt + 
## mod.tree.shannon3:     (1 | quad.unique)
## mod.tree.shannon1: shannon ~ site + census + pesticide * exclosure + (1 | quad.unique)
##                   Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.tree.shannon2 12 516.04 573.03 -246.02   492.04              
## mod.tree.shannon3 13 517.96 579.69 -245.98   491.96 0.0892      1
## mod.tree.shannon1 15 521.59 592.83 -245.80   491.59 0.3612      2
##                   Pr(>Chisq)
## mod.tree.shannon2           
## mod.tree.shannon3     0.7652
## mod.tree.shannon1     0.8348
plot(mod.tree.shannon2, type=c('p', 'smooth'))

qqPlot(resid(mod.tree.shannon2))

sjp.lmer(mod.tree.shannon2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.tree.shannon2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.tree.shannon2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
##    Data: subset(divers_tspest, growth.form == "tree")
## 
## REML criterion at convergence: 543.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4107 -0.5661  0.0429  0.6062  2.9724 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.06964  0.2639  
##  Residual                0.06892  0.2625  
## Number of obs: 853, groups:  quad.unique, 239
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.81838    0.04553  17.975
## site1        0.09859    0.02740   3.598
## site2        0.03608    0.02724   1.324
## census16sp  -0.69200    0.02888 -23.958
## census17fa   0.23466    0.02410   9.735
## census17sp  -0.10985    0.02412  -4.555
## pesticideF   0.12179    0.05440   2.239
## pesticideI   0.03059    0.05487   0.558
## pesticideC  -0.06755    0.05471  -1.235
## exclosure1  -0.07763    0.03870  -2.006
## Invsimpson
divers_tspest1 <- subset(divers_tspest, simpson != Inf)
mod.tree.simpson <- lmer(simpson ~ site + census + pesticide + exclosure + (1|quad.unique), 
                         data=subset(divers_tspest1, growth.form=='tree'))

plot(mod.tree.simpson, type=c('p', 'smooth'))

qqPlot(resid(mod.tree.simpson))

sjp.lmer(mod.tree.simpson, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.tree.simpson, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.tree.simpson, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## simpson ~ site + census + pesticide + exclosure + (1 | quad.unique)
##    Data: subset(divers_tspest1, growth.form == "tree")
## 
## REML criterion at convergence: 1764.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8480 -0.5775 -0.0833  0.5437  3.5601 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.3053   0.5525  
##  Residual                0.3258   0.5707  
## Number of obs: 815, groups:  quad.unique, 238
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.16286    0.09681  22.340
## site1        0.23976    0.05823   4.118
## site2        0.02483    0.05789   0.429
## census16sp  -0.99589    0.06912 -14.407
## census17fa   0.59961    0.05263  11.394
## census17sp  -0.14007    0.05265  -2.660
## pesticideF   0.24869    0.11558   2.152
## pesticideI   0.12667    0.11689   1.084
## pesticideC  -0.06748    0.11609  -0.581
## exclosure1  -0.21902    0.08230  -2.661
#### shrub
# shannon diversity
mod.shrub.shannon1 <- lmer(shannon ~ site + census + pesticide*exclosure + (1|quad.unique), 
                           data=subset(divers_tspest, growth.form=='shrub'))
mod.shrub.shannon2 <- lmer(shannon ~ site + census + pesticide + exclosure + (1|quad.unique),
                           data=subset(divers_tspest, growth.form=='shrub'))
mod.shrub.shannon3 <- lmer(shannon ~ site + census + pesticide + exclosure + A.sum_curt +
                             (1|quad.unique), data=subset(divers_tspest, growth.form=='shrub'))
anova(mod.shrub.shannon1, mod.shrub.shannon2, mod.shrub.shannon3)
## refitting model(s) with ML (instead of REML)
## Data: subset(divers_tspest, growth.form == "shrub")
## Models:
## mod.shrub.shannon2: shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
## mod.shrub.shannon3: shannon ~ site + census + pesticide + exclosure + A.sum_curt + 
## mod.shrub.shannon3:     (1 | quad.unique)
## mod.shrub.shannon1: shannon ~ site + census + pesticide * exclosure + (1 | quad.unique)
##                    Df     AIC     BIC logLik deviance  Chisq Chi Df
## mod.shrub.shannon2 12 -191.34 -135.85 107.67  -215.34              
## mod.shrub.shannon3 13 -191.40 -131.29 108.70  -217.40 2.0587      1
## mod.shrub.shannon1 15 -186.12 -116.76 108.06  -216.12 0.0000      2
##                    Pr(>Chisq)
## mod.shrub.shannon2           
## mod.shrub.shannon3     0.1513
## mod.shrub.shannon1     1.0000
plot(mod.shrub.shannon2, type=c('p', 'smooth'))

qqPlot(resid(mod.shrub.shannon2))

sjp.lmer(mod.shrub.shannon2, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.shrub.shannon2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.shrub.shannon2, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## shannon ~ site + census + pesticide + exclosure + (1 | quad.unique)
##    Data: subset(divers_tspest, growth.form == "shrub")
## 
## REML criterion at convergence: -167.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4242 -0.2175 -0.0030  0.1953  3.9673 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.22503  0.4744  
##  Residual                0.01582  0.1258  
## Number of obs: 753, groups:  quad.unique, 195
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  0.6630352  0.0778860   8.513
## site1        0.1674967  0.0470239   3.562
## site2        0.0603374  0.0484612   1.245
## census16sp  -0.0247005  0.0130640  -1.891
## census17fa   0.0286765  0.0130150   2.203
## census17sp   0.0002409  0.0129706   0.019
## pesticideF  -0.0187711  0.0984434  -0.191
## pesticideI  -0.0345274  0.0973471  -0.355
## pesticideC  -0.0640912  0.0963271  -0.665
## exclosure1  -0.0697174  0.0688905  -1.012
## Invsimpson
mod.shrub.simpson <- lmer(simpson ~ site + census + pesticide + exclosure + (1|quad.unique), 
                          data=subset(divers_tspest1, growth.form=='shrub'))

plot(mod.shrub.simpson, type=c('p', 'smooth'))

qqPlot(resid(mod.shrub.simpson))

sjp.lmer(mod.shrub.simpson, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.shrub.simpson, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.shrub.simpson, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## simpson ~ site + census + pesticide + exclosure + (1 | quad.unique)
##    Data: subset(divers_tspest1, growth.form == "shrub")
## 
## REML criterion at convergence: 867.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7200 -0.2548 -0.0147  0.1818  3.7375 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.88945  0.9431  
##  Residual                0.06534  0.2556  
## Number of obs: 747, groups:  quad.unique, 192
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  2.114774   0.155246  13.622
## site1        0.320558   0.094192   3.403
## site2        0.053499   0.097229   0.550
## census16sp  -0.051347   0.026575  -1.932
## census17fa   0.061410   0.026476   2.319
## census17sp  -0.002876   0.026533  -0.108
## pesticideF  -0.066060   0.198170  -0.333
## pesticideI  -0.054336   0.194706  -0.279
## pesticideC  -0.121789   0.191635  -0.636
## exclosure1  -0.092261   0.138194  -0.668

When trees and shrubs are pooled in the diversity test, pesticide had no effect on diversity. Exclosure lower diversity. When trees and shrubs are estimated seperately, fungicide higher diversity of tree seedlings.